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SECTION  I 
INTRODUCTION 


For  aircraft  weapon  systems  which  deliver  unguided  armament  (guns, 
rockets,  bombs),  small  improvements  in  delivery  accuracy  can  result  in 
substantial  reductions  of  weapon  dispersion  [  1] .  The  net  result  is  a  reduced 
number  of  sorties  required  per  target  and  hence  less  exposure  of  aircraft  to 
enemy  defenses.  Consequently,  the  analysis  of  precision  weapon  delivery 
has  become  an  area  of  extreme  interest.  Two  approaches  to  the  problem 
have  evolved.  The  first  is  to  analyze  the  precision  weapon  delivery  task  by 
detailed  simulation  of  the  tactical  situation  and  the  second  is  to  investigate 
the  effects  of  initial  condition  errors,  release  point  errors,  and  gust  dis¬ 
turbance  errors  by  propagating  them  along  the  weapons  trajectories  to  im¬ 
pact.  For  the  latter  approach,  a  method  of  analysis  is  needed  which  con¬ 
siders  the  effects  of  flight  control  parameters,  airframe  dynamics,  measure¬ 
ment  errors  and  gust  disturbances  on  the  release  parameters  and  then 
propagates  these  release  errors  to  impact  to  obtain  a  meaningful  measure 
of  weapon  system  performance. 

Such  an  analysis  tool  can  then  be  used  in  connection  with  the  former 
approach  to  establish  critical  measurement  requirements  as  well  as  to 
evaluate  the  effects  of  various  control  points  and  methods  prior  to  large- 
scale  simulation  or  flight  tests. 

In  this  study,  a  dynamic  precision  weapon  delivery  system  model,  a 
method  of  performance  analysis,  and  a  set  of  computer  programs  are  devel¬ 
oped  to  evaluate  effects  of  various  process  parameters  on  overall  weapon 
system  performance,  with  particular  emphasis  on  flight  control  parameters, 
airframe  dynamics,  measurement  points  and  errors,  and  gust  environments. 
The  model  of  the  process  assumes  precision  control  without  considering  the 
dynamics  of  the  operator  as  a  control  element.  The  controller  is  in  the  form 
of  an  optimal  controller  consisting  of  an  optimal  estimator  and  optimal  feed¬ 
back  gains.  The  model  is  general  enough,  however,  to  consider  existing 
controllers  as  well  as  the  operator  in  the  loop.  For  the  latter  case  (i.  e. , 
pilot  weapon  delivery)  mathematical  models  of  the  operator  have  to  be  in¬ 
corporated  into  the  present  model  [2,  3,  4,  5],  Since  the  operator  tracking 
error  is  one  of  the  major  contributors  to  impact  dispersion,  the  piloted 
delivery  deserves  separate  attention.  The  model  is  also  flexible  enough  for 
considering  alternate  airframe  dynamics /control  points /measurement  system 
combinations.  The  developed  armament  model  is  for  the  delivery  of  iron 
bombs.  It  is  sufficiently  general  to  treat  a  variety  of  dive  angles  and  release 
altitudes  and  bomb  characteristics.  The  delivery  of  other  weapons  can  be 
considered  also  by  minor  modifications. 

The  impact  covariance  and  the  circular  error  probable  (CEP)  are  used 
as  measures  of  delivery  performance. 
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A  simplified  variance  analysis  is  given  in  [5]  by  ignoring  cross- 
covariance  terms.  The  elimination  of  the  cross -variance  terms  was  pur¬ 
sued  in  [2]  by  choosing  random  variables  which  are  uncorrelated. 

In  this  work  the  full  impact  covariance  matrix  is  developed  and  retained 
in  the  analysis  using  the  complete  dynamical  model.  The  concepts  of  half¬ 
probability  circle  and  circular  error  probable  are  extended  to  half-probability 
ball  and  spherical  error  probable.  These  are  thought  to  be  more  meaningful 
performance  measures  in  air-to-air  delivery. 

The  complete  report  is  divided  into  three  volumes.  Volume  I  contains  the 
works  on  the  weapon  delivery  system  modeling  and  optimization.  Volume  II 
documents  the  programs  which  implement  the  analysis  developed  in  Volume  I. 
Volume  III  contains  a  demonstration  example  to  illustrate  how  these  programs 
are  used. 

In  this  volume  the  presentation  begins  with  a  brief  description  of  the  wea¬ 
pon  delivery  process  and  the  approach  taken  for  the  performance  analysis  of 
the  overall  system.  In  addition,  the  synthesis  of  the  optimal  weapon  delivery 
controller  is  outlined  and  the  overall  organization  of  the  analysis  program  is 
given.  In  Section  III  the  development  of  a  mathematical  model  for  the  six- 
degree-of-freedom  motion  of  aircraft  and  weapon  is  presented.  In  Section  IV 
the  development  of  the  total  force  and  moment  system  for  the  aircraft  and 
weapon  moving  in  an  unsteady  airmass  is  given.  The  wind  model  is  developed 
in  that  section  also.  The  development  of  the  measurement  system  model  is 
presented  in  Section  V.  The  nonlinear  observation  geometry  as  well  as  sensor 
dynamics  are  considered.  In  Section  VI  the  process  of  linearization  is  treated. 
The  transformation  of  the  perturbation  states  is  presented  in  that  section  also. 
The  development  of  the  nominal  states  and  parameters  (trimming)  for  the 
linearization  is  discussed  in  Section  VII.  The  methods  of  algebraic  as  well 
as  the  autopilot  trim  are  presented.  The  performance  measure  development 
for  the  analysis  and  design  of  the  weapon  delivery  controller  is  given  in 
Section  VIII.  The  circular  error  probable  performance  index  is  generalized 
to  the  spherical  error  probable  for  the  air-to-air  weapon  delivery.  In  Sec¬ 
tion  IX  a  method  for  nonstationary  optimal  weapon  delivery  controller  design 
is  presented.  Both  deterministic  and  stochastic  disturbances  are  considered. 
This  is  followed  by  the  presentation  of  a  stationary  controller  design  method 
in  Section  X.  The  iterative  solutions  of  the  Lyapunov  as  well  as  the  Riccati 
equations  are  given  in  that  section  also.  Section  XI  summarizes  the  analysis 
and  modeling  work  and  lists  recommendations  for  additional  areas  of  study 
and  extensions. 

The  reduced  controller  (fixed-form)  design  method  is  given  in  Appen¬ 
dix  I.  The  development  of  the  fire  control  equations  for  the  nominal  release 
time  is  given  in  Appendix  II. 
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SECTION  II 

PERFORMANCE  ANALYSIS  AND  OPTIMAL  DESIGN 
OF  WEAPON  DELIVERY  PROCESSES 

This  section  provides  an  overview  of  the  approach  taken  for  the  system 
analysis  and  the  optimal  design  of  weapon  delivery  processes.  First,  a 
brief  description  of  the  air-to-ground  weapon  delivery  processes  is  pre¬ 
sented;  then  the  mathematical  models  of  the  subsystems  and  the  overall 
delivery  system  are  described.  Subsequently  the  optimal  perturbation 
control  of  the  delivery  model  is  discussed,  and  finally,  the  overall  organi¬ 
zation  of  the  "Armament  Delivery  Analysis  Programming  System  (ADAPS)" 
is  given. 


AIR-TO-GROUND  DELIVERY  PROCESSES 

Air-to-ground  delivery  is  currently  accomplished  with  a  limited  number 
of  well-defined  attack  maneuvers,  each  designed  for  a  particular  tactical  situ¬ 
ation  (i.  e. ,  for  given  target  type  and  defenses,  aircraft  type,  armaments  and 
electronic  aids),  figure  1  shows  some  of  the  basic  dive-bombing  trajectories 
used  mostly  for  manual  (iron-sight)  delivery  [7],  In  this  delivery  technique, 
the  pilot  flies  a  fairly  consistent  preplanned  flight  path  which  brings  the  air¬ 
craft  to  an  initially  set  release  condition,  e.g. ,  release  altitude,  ground 
speed  and  dive  angle.  The  weapon  is  released  when  these  three  conditions 
are  satisfied  simultaneously  in  the  ideal  cases  and  the  target  passes  through 
the  center  of  a  depressed  reticle  sight.  The  depressed  reticle  sight  display 
system  is  used  most  often  for  target  tracking  during  conventional  air-to- 
ground  weapon  delivery. 

The  aim  dot  (pipper)  of  the  reticle  image,  which  the  pilot  attempts  to 
position  onto  the  target  by  steering  the  aircraft,  is  "depressed"  below  the 
velocity  vector  of  the  aircraft  as  shown  in  Figure  2  in  order  to  display  the 
nominal  bomb  impact  point  corresponding  to  the  preselected  nominal  release 
condition. 

The  setting  of  the  sight  angle  is  determined  from  the  nominal  trajectory 
line-of-sight  angle,  yios.  The  sight  depression  angle,  T),  is  the  angle  of 
the  line  of  sight  below  the  fuselage  line  (x-axis)  of  the  aircraft.  The  pilot 
determines  from  tabulated  data  the  correct  fixed  sight  depression  angle  for 
the  selected  pickle  (release)  conditions,  and  adjusts  the  sight  depression 
mechanism  to  this  setting.  He  dives  towards  the  target  in  such  a  way  as  to 
achieve  the  nominal  release  conditions  at  the  time  the  pipper  is  on  the  target. 

When  using  fire  control  computers,  achieving  the  same  accurate  dive 
conditions  is  not  necessary.  Instead,  the  target  is  held  on  the  center  of  an 
"undepressed"  reticle  sight  throughout  the  dive,  and  the  weapon  system  is 
pickled  at  the  desired  release  altitude.  This  initiates  the  weapon  system  com¬ 
puter  for  automatic  release. 
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The  automatic  release  capability  increases  the  number  of  possible  de¬ 
livery  maneuvers.  Figure  3  shows  a  typical  example,  the  "high- approach 
dive-toss"  maneuver  with  automatic  release  [  8] .  After  target  acquisition, 
this  maneuver  may  be  broken  into  four  phases: 

•  Dive  --  The  aircraft  is  aligned  with  ihe  target  and  flown 

to  a  desired  release  altitude,  ground  speed,  and  dive  angle. 

The  automatic  weapon  release  computer  set  (WRCS)  is 
(pickled  at  the  desired  altitude)  for  starting  the  pull-up. 

This  initiates  the  weapons  system  computer  for  auto¬ 
matic  release. 

•  Pull-up  —  Following  the  release  or  pickle  command,  a 
constant-g  pull-up  is  initiated.  This  serves  the  obvious 
purpose  of  reducing  exposure  to  enemy  defenses  and  also 
provides  for  clean  aircraft-bomb  separation.  With  fire 
control,  weapon  release  will  occur  automatically  during  the 
pull-up  maneuver. 

•  Release  Transient  —  Weapons  are  commonly  released  from 
their  mounting  racks  with  explosive  charges  that  inject  them 
into  die  airflow  near  wings  or  fuselage  with  uncertain  linear 
and  angular  velocities.  The  resulting  short-lived  but  poorly 
understood  transients  establish  initial  conditions  for  the 
weapon's  free  fall. 

•  Weapon  Free-Fall  —  The  weapon  follows  a  ballistic  trajectory 
toward  the  target.  This  may  be  single-stage,  as  for  iron 
bombs,  or  multi-stage,  as  for  dispenser-type  weapons.  This 
work  deals  primarily  with  free-falling  weapons,  subject  only 
to  aerodynamic  and  gravity  forces.  The  analysis  program, 
however,  is  sufficiently  general  to  accept  nonballistic  tra¬ 
jectories  for,  say,  guided  bombs  or  missiles. 

While  the  dive-toss  is  a  specialized  attack  maneuver,  it  can  easily  be 
parameterized  with  respect  to  acquisition  altitude,  dive  angle,  speed, 
release  altitude  and  pull-up  g's  to  generate  various  other  maneuvers.  For 
example,  letting  the  dive  angle  y  vanish  and  executing  a  zero-g  increment 
pull-up  (and  depressing  the  sight  reticle  appropriately)  leads  to  the  standard 
"low-level  approach  and  laydown"  maneuver.  Similarly,  y  =  0  with  positive-g 
pull-up  and  automatic  weapon  release  gives  a  "loft  bombing"  maneuver. 

Because  of  this  inherent  generality,  the  discussions  to  follow  are  built 
around  the  dive-toss  maneuver  as  a  typical  nominal  trajectory. 
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Figure  3.  High- Approach  Dive  Toss  (with  automatic 
release) 
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MATHEMATICAL  MODELS 


The  weapon  delivery  processes  in  this  work  are  treated  as  nonlinear 
stochastic  phenomena.  This  is  both  realistic  and  computationally  attractive. 
For  each  phase  of  the  attack  maneuver,  the  process  is  linearized.  CEP  per¬ 
formance  is  evaluated  with  standard  covariance  analyses,  and  optimization 
is  accomplished  with  linear-quadratic  control  theory.  No  repeated  runs  or 
Monte  Carlo  methods  are  required.  Discussed  below  are  the  models  appli¬ 
cable  to  each  trajectory  phase  --  dive,  pull-up,  release  transient,  and  weapon 
free-fall. 


Dive  and  Pull-up  Phases 


The  mathematical  model  for  the  dive  and  pull-up  phases  includes  aircraft 
equations  of  motion,  measurement  equations  with  measurement  noises  and 
biases,  and  wind  model  equations  for  gusts  and  mean  winds.  The  model  also 
includes  simple  fire  control  equations  which  define  the  nominal  pull-up  and 
release  times.  The  overall  system  contains  both  linear  and  nonlinear  dyna¬ 
mics  (Figure  4). 


Airframe  Model  --To  handle  arbitrary  trajectories,  various  control  points, 
and  methods,  the  aircraft  mode1  is  described  by  a  set  of  nonlinear,  time- 
varying  differential  equations  of  the  form 

x  =  f(x  ,  x  ,  yT,  y.,  w,  v) 

where 


state  vector  of  the  aircraft 

vector  of  effective  thrust  inputs  to  the  aircraft 

vector  of  effective  surface  deflection  inputs  to  the  aircraft 

vector  of  disturbances  (gusts)  on  the  aircraft 

vector  of  deterministic  inputs  (mean  winds)  to  the 
aircraft 


Actuator  and  Thrust  System  Model  —  To  allow  variations  in  the  effectiveness 
of  thrust  and  aerodynamic  surfaces,  the  actuator  and  thrust  system  outputs 
are  assumed  to  be  nonlinear  functions  of  their  states. 


The  actuator  and  thrust  dynamics  are  modeled  as  follows: 


x.  =  +  G.u.  +  K  ,y „ 

6  6  6  6  6  y  5  m 

Xrp  ~  FfJ*Xrp  4"  GrpUq,  +  Ky  rp  ^ 


with  nonlinear  outputs 

y8  =  W 

yT  *  hr<JtT> 


where 


=  actuator  system  state  vector 


x, p  =  thrust  system  state  vector 

y.  =  effective  aerodynamic  surface  deflection  output  vector 
o 

y  =  effective  thrust  magnitude  and  deflection  output  vector 

K  K  T  are  arbitrary  feedback  gains  and  y  is  the  measured 
y  1  signals  as  modeled  below. 

Measurement  System  Model  --To  allow  various  measurement  points  and 
methods ,  the  measurement  system  model  is  represented  by  a  set  of  linear 
differential  equations  (i  e. ,  sensor  dynamics)  with  linear  and  nonlinear  in¬ 
puts.  The  nonlinear  inputs  correspond  to  the  observation  geometry  of  the 
state  measurements 

Xm  ~  Fm  xm+  ^mp  ypm  +  ^mp  ^pm+^mT  yTm  +  ^mfi 
+  Gmnm 


y  =  M  x  +  D _ y  +  n 

J  m  mm  me  J  e  m  'a 


where 


xm  =  measurement  system  state  vector 

Y  =  sensed  signal  vector  for  aircraft  states 
J  pm  6 

• 

y  sensed  signal  vector  for  aircraft  state  derivatives 

y  _  =  sensed  signal  vector  for  the  thrust  magnitude  and  thrust 
m  deflection  states 

y^m=  sensed  signal  vector  for  the  actuator  states 

^m  =  w^te  no*se  inPut  vector  he  measurement  system 

ym  =  measured  signals 


0 


ye  =  output  vector  from  the  bias  error  process 
T)  =  additive,  white  measurement  noise 

St 

The  sensed  signal  vectors  for  the  aircraft  states  are  assumed  to  be  non¬ 
linear  (i.e. ,  observation  geometry)  functions  of  the  states  and  its  derivatives: 
♦  •  • 

ypm  =  h'V  y 

fpm  ■  h<y 

The  sensed  signal  vectors  for  the  thrust  and  actuator  states  are  assumed 
to  be  linear  functions  of  the  thrust  and  actuator  states 

y  _  =  H„  x  „ 

J  Tm  Tm  T 


=  H5mx6 

The  measurement  system  model  not  only,  contains  those  measurements 
for  controlling  the  aircraft  in  flight  (such  as  gyro  and  accelerometer  outputs 
used  for  augmenting  pilot  control  and  for  flexure  control),  but  also  those  mea¬ 
surements  necessary  to  estimate  the  states  for  the  fire  control  equations 
(such  as  radar  range,  azimuth  and  elevation,  and  altimeter  outputs  to  deter¬ 
mine  time  of  weapons  release).  Error  sources  in  the  measurement  system 
are  also  included.  Examples  of  these  are  gyro  biases,  accelerometer  noise 
due  to  engine  vibrations,  and  radar  angle  errors  due  to  misalignment. 

Linear  Equations  --For  small  perturbation  analysis  the  nonlinear  part  of  the 
overall  system  (i.  e. ,  the  nonlinear  equations  of  motion,  the  effective  input 
equations  and  the  observation  equations)  is  linearized  about  the  nominal 
flight  path  to  the  form 


•  bf 
gx  =  ~ 
P  3xr 


gX  +  vT- 

P  dxr 


6X  +|3L 

p  |9yT 


3hT  * 

5xt 


6yx  +{^)6yB + 


pm  dx 


gx  + 
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The  measurement  perturbation  equations  become 


fix  =  F  fix  +  G.  5y_  +  G  fiy  +G  mfiyrp^, 

m  mm  tup  ym  mp  J  pm  mi  Tm 


+  Gmfi  6yfim  +  Gm  11  m 
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the  state  perturbation  equations  become 

%  ’  Fp  6  V  V  V  V  HTp5xT  +  Gp6  H6p  6lte  +  GpwW  +  Gpv'' 


fix  =  F  fix  +  (G  H  fix  +  G  H  )  fix  +  G  H  fix 
m  m  m  mp  pm  p  mp  pm7  p  mp  pm  p 


+  GmT  HTm  6xT  +  Gmfi  Hfim  6xfi  +  Gm  ^  m 


The  wind  gust  model  with  the  Dryden  spectrum  is  described  by 
x 


w 


F  (t)  x  +  G  n 
w  w  w  w 


w  =  H  x 
wp  w 


where 


xw  =  gust  filter  state 


w  =  gust  filter  output 

=  white  noise  input  to  the  filter 


The  measurement  bias  error  vector  xe(t)  is  modeled  as  a  differential 
equation  with  a  long  time  constant  so  that  it  remains  approximately  constant 
over  the  time  of  the  weapon  delivery  maneuver.  That  is 
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F  x  +  G  ri 
e  e  e  e 


H  x 
em  e 


bias  error  vector  state 

error  system  output 

white  noise  input  to  error  filter 


The  matrix  Ge  is  selected  so  that  the  steady-state  rms  values  of  xg  are 
equal  to  the  rms  biases. 


Defining 


x 

y 

u 

r\ 


col  [*m,  *5,  *T.  Xp,  Xw.  xel 


Jm 

=  col  [Ug,  Uj.3 

■  “l  C  V  V  V  V 


the  augmented  mathematical  model  is  obtained  in  the  form  of 


Cl-F(t)]  k  *  F(t)x  +  Gu(t)u  +  G-  (t)v  +  G^ 


y  =  M(t)x  +  Dr| 


where  the  matrices  F,  Gu,  G-,  G^,  M  and  D  are  given  as  follows: 


F  = 


m 


GmeHem 


G  _H_  G  H  +  G  H  0 
mT  Tm  mp  pm  mp  pm 
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This  finishes  the  dive  phase  of  the  modeling.  Figure  4  shows  the  overall 
system  linearized  state  diagram. 

Release  Phase 

For  the  automatic  release,  equations  predicting  the  nominal  pull-up 
time  and  the  release  time  are  needed.  They  are  in  the  form  of 


13 


where 


tp  =htp[x(t)j,  t<tp 
tr  =htr[x(tp).  x(t)],  tp<t<tr 


t  =  nominal  pull-up  time 
P 

tr  =  nominal  time  of  weapon  release 
x  =  state  vector  of  the  nonlinear  model 

These  are  linearized  at  the  nominal  release  time  t  =  tp  yielding 


6x(tp) 


SV 

6tr =  My  C“(tp)]  + 

where  prime  indicates  the  transpose. 


ah. 


dx 


|  [6x(tr>] 


If  the  estimate  of  perturbation  state  6x  is  available  instead  of  6x  itself 
then  in  the  above  equations  6x  is  replaced  by  6x.  These  equations  define  the 
timing  errors  in  the  release  time.  Computation  delay  and  rack  relay  delays 
contribute  additional  timing  errors.  These  are  modeled  as  independent  addi¬ 
tive  random  variables. 


Release  Transient  Phase 


The  release  transient  phase  is  defined  here  as  the  period  from  the  actual 
release  time  of  bomb  to  the  time  when  it  leaves  the  region  of  the  wings.  Al¬ 
though  this  phase  is  not  modeled  in  detail  in  this  work,  it  has  an  important 
influence  on  the  weapon  delivery  performance.  The  records  in  some  cases 
show  that  when  the  bomb  is  ejected  from  the  wing  rack,  the  wing  moves  up 
(flexure)  while  the  bomb  follows  a  trajectory  without  safe  separation.  Because 
of  the  inability  to  predict  bomb  aerodynamics  and  stability  in  the  region  of 
the  wing,  experimental  data  are  used  for  the  description  of  the  release 
transient  phase.  At  the  present  time,  the  U.S.  Air  Force  has  a  program 
called  "SEEK  EAGLE"  [9, 10],  which  is  both  analytically  and  experimentally 
studying  the  release  problem.  In  this  work,  the  release  transient  phase  is 
taken  into  account  by  introducing  at  release  an  independent,  additive,  stochastic, 
initial-condition  error. 
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The  state  of  bomb  at  release  is  described  by 

xb  ■  hb  <y 

The  perturbation  state  of  bomb  at  nominal  release  is  then  given  by 


where 


•W  =  Hb  6xp  (*r> 


Hb =  ar 


The  state  of  bomb  immediately  following  release  transient  is  given  by 


«V*r+>  =  Hb  [6xp(lr)  +  fr  6tr]+  Hr  ?, 


where 


fr  =  xp(tr} 

5t  =  release  time  error 
r 

Hp  =  release  transient  error  input  matrix 

£  =  release  transient  error 

3r 

Weapon  Free-Fall  Phase 

In  evaluating  weapon  delivery  performance,  translating  the  dispersion 
errors  at  bomb  release  to  the  impact  of  the  bomb  on  the  target  and  introducing 
errors  which  occur  during  the  bomb  trajectory  is  necessary.  Two  possibili¬ 
ties  exist  for  satisfying  this  requirement  —  bomb  tables  or  trajectory  compu¬ 
tation.  As  trajectory  computation  allows  treatment  of  various  weapons  and 
permits  introduction  of  disturbances  such  as  winds  and  uncertain  bomb  param¬ 
eters  during  the  trajectory,  this  method  has  been  chosen  in  this  work. 

The  bomb  model  is  described  by 


xfa  =  f(xb,  xb,  u,  w,  v) 


where 
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x,  =  state  vector  of  the  bomb 
b 

u  =  input  vector  due  to  bomb  imperfections 

w  =  vector  of  disturbances  (gusts)  on  the  bomb 

v  =  vector  of  deterministic  inputs  (mean  winds)  to  the  bomb 

For  small-perturbation  analysis,  the  nonlinear  equations  are  linearized 
about  the  nominal  free-fall  flight  path.  This  yields 

6xb  =  Fb(t)  6xb  +  Fb(t)  6kb  +  Gbw(t)  w  +  Gbyv  +  Guu 

The  perturbation  state  of  bomb  at  the  impact  on  the  horizontal  plane  is 
given  by 

6x(tf )  =  6x  (tf)  +  fb(tf)  6tf 


where 

6x(tj)  =  the  perturbation  state  of  bomb  at  the  nominal  impact  time  t^ 
W  =  xb(tf) 

and 

6h(t.) 

St  =  - -I— 

0tf  H(tf) 

h,  being  the  altitude  state  component  of  the  bomb. 

After  having  modeled  the  delivery  process  from  the  target  acquisition  to 
the  impact,  various  performance  measures  can  be  defined  for  the  analysis 
and  design  of  the  delivery  system  as  presented  in  the  following  subsection. 


OPTIMAL  CONTROL  OF  PERTURBATION  MODELS 

The  mathematical  models  just  discussed  provide  small  perturbation 
descriptions  of  the  weapon  delivery  process.  They  are  incomplete,  however, 
in  that  the  control  variables  ug  and  u-p  for  the  dive  and  pull-up  phases  of  the 
trajectory  are  undefined.  Adding  arbitrary  (linear)  controllers  to  the  mathe¬ 
matical  model  is  a  simple  matter,  and  the  developed  aircraft  model  provided 
this  option.  Its  utility  lies  in  quick  performance  evaluations  of  all  kinds  of 
specific  bomb  systems  and  control  schemes. 

Tools  are  needed  for  evaluating  intrinsic  properties  of  various  constraint 
configurations  to  answer  questions  such  as; 
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•  For  a  fixed  complement  of  sensors  and  a  given  set  of  control 
inputs  (control  points),  what  minimum  CEP  is  attainable  within 
these  constraints? 

•  What  maximum  CEP  improvements  does  adding  a  single  sensor 
to  the  above  configuration  give? 

•  What  performance  penalty  exists  for  replacing  a  sensor  with  a 
noisier,  but  cheaper,  version? 

These  questions  call  for  performance  analyses  of  optimal  controllers. 

A  large  number  of  possible  optimal  control  problems  can  be  formulated 
for  the  weapon  delivery  process  —  for  example,  computing  optimal  attack 
trajectories.  This  has  already  been  avoided  by  saying  that  these  are  largely 
predetermined  by  the  tactical  situation.  Another  problem  is  that  of  optimally 
controlling  velocity  deviations  and  target  deviations  from  the  reticle  of  the 
bomb  sight  during  the  dive  phase  of  attack,  followed  either  by  an  open-loop 
or  optimally  controlled  pull-up.  Assuming  the  steady-state  operation,  the 
performance  index  for  this  formulation  takes  the  form 

h  -  E  {qn  6V2  +  q22  %2  +  q33  ey2  +  u‘ '  Ru} 
where 

6V  =  velocity  deviation  from  nominal 

=  vertical  deviation  of  target  from  center  of  reticle 

e  =  lateral  deviation  of  target  from  center  of  reticle 

i# 

q..,  q22»  q33,  R  -  weighting  coefficients  and  prime  indicates  the 

transpose 

The  solution  mimics  a  pilot's  own  efforts  to  hold  target  alignment  and 
velocity  during  his  dive.  The  weights  q^,  R  can  be  chosen  via  the  iterative 
methods  of  quadratic  equivalence  [111 

Still  another  optimization  problem  which  is  followed  in  this  work  involves 
direct  minimization  of  the  HPA.  This  can  be  done  with  a  performance  index 
obtained  in  the  following  manner. 

Strictly  speaking,  HPA  is  the  area  of  a  0.  5  probability  circle  centered 
at  the  mean  impact  point;  CEP  is  its  radius.  For  normal  distributions  with 
small  cross  correlations,  this  area  can  be  closely  approximated  by 

HPA  =  [qx(7x2  <tf)  +  qyoy2  <tf)] 

where  Ox  and  Oy  are  downrange  and  crossrange  standard  deviations,  and  qx, 
qy  are  weightings  which  depend  on  the  ratio  Ox/ ay. 


J?FZW?x f^jw XtansKRZEim?? r^’-'n  -• 
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This  expression,  in  turn,  can  be  written  in  terms  of  the  bomb  release 
covariance  matrix  X(tr+)  using  the  perturbance  equations  of  the  bomb 

HP  A  =  tr  {X(tr+)  Qr  +  XfQf } 


where 

«. 

Q, 


tr  = 


propagation  weighting  matrix  for  release  errors 
propagation  weighting  matrix  for  the  forced  response 
impact  covariance  response  due  to  wind  gusts, 
trace  operator 


The  expression  given  above  shows  that,  with  this  performance  measure, 
performance  analysis  of  a  weapon  delivery  process  reduced  to  standard  linear 
covariance  analysis. 

To  directly  minimize  HP  A,  therefore,  the  performance  index  should  be 

t 

J0  =  HP  A  +  rr  tr  fu'  R  u}dt 


( 


Note  that  this  is  a  performance  index  with  the  terminal  cost,  penalizing 
errors  at  the  nominal  release  time  tr-  This  means  that  the  optimal  controller 
will  be  time  varying  even  if  the  system  dynamics  are  stationary.  This  makes 
for  expensive  analysis.  If  stationary  dynamics  are  adequate,  modifying  the 
performance  index  such  that  it  yields  a  stationary  controller  is  desirable. 

A  steady-state  version  of  Jo  does  just  that: 

J3  =  HP  A  +  tr  {R  U} 


where 

U  =  E  {u  u'} 

An  important  distinction  exists  between  the  optimization  problems  based 
on  J2  and  Jo  and  the  optimization  problem  based  on  the  earlier  index  J,. 
Solutions  01  the  Jj  problem  depend  on  the  type  of  sight  in  the  aircraft  — 
whether  it  is  fixed,  or  drift  stabilized,  or  pitch  stabilized,  etc.  Solutions  of 
the  J2,  J3  problems,  on  the  other  hand,  completely  bypass  the  sighting  sys¬ 
tem.  They  depend  only  on  the  basic  constraint  configuration,  i.  e. ,  on  the 
signals  available  for  measurement  and  on  the  control  points  available  for 
manipulation. 
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OVERALL  ORGANIZATION  OF  ADAPS 


The  various  subroutines  implementing  the  above  model  provide  the  capa¬ 
city  to  analyze  weapon  delivery  as  a  general  linear  time-varying,  stochastic 
process  or  to  analyze  it  as  a  much  simplified  process  that  is  stationary  during 
each  of  its  phases.  The  one  extreme  offers  fidelity  to  the  physical  situation, 
while  the  other  offers  low  computing  costs  and  the  possibility  of  many  analysis 
iterations.  By  using  the  program  organization  shown  in  Figure  5,  both  ex¬ 
tremes  (as  well  as  the  many  possibilities  in  between)  are  readily  attainable. 

In  this  organization,  the  individual  subroutines  are  accessible  from  a  main 
program  with  which  they  share  common  memory.  They  communicate  with 
each  other  within  the  groups  indicated.  Optional  inputs  are  provided  to  cover 
the  various  special  possibilities  discussed  above.  A  detailed  description  of 
the  overall  organization  of  ADAPS  is  presented  in  the  Section  II  of  Volume  II. 
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A  typical  analysis  proceeds  as  follows: 

•  I  -  Linear  Data  Generation 

1.  Read  input  data  for  attack  maneuver,  read  nonlinear  air¬ 
craft  aerodynamics . 

2.  Trim  aircraft,  and  fly  it  to  obtain  nominal  trajectory  up  to 
nominal  release  altitude. 

3.  Linearize  the  aircraft  equations  of  motion  numerically  at 
specified  time  points  during  flight.  Write  on  tape. 

4.  Read  input  data  for  nonlinear  weapon  aerodynamics. 

5.  Using  the  release  conditions,  generate  free-fall  trajectory 
by  the  nonlinear  weapon  model. 

6.  Linearize  the  weapon  equations  of  motion  numerically  at 
specified  time  points  along  the  free-fall  trajectory,  write 
on  tape. 

•  II  -  Optimization 

1.  Generate  the  propagation  weighting  matrix  using  linear 
weapon  data. 

2.  Input  control  points,  measurement  points,  measurement 
variances. 

3.  Choose  the  type  of  optimization,  and  obtain  controller  gains, 
estimator  gains,  total  system  covariance  at  release. 

•  in  -  Performance  Evaluation 

1.  Propagate  the  release  covariance  to  impact  using  performance 
evaluator. 

2.  Compute  impact  covariance  matrix,  CEP  performance 
measure  and  the  variance  contribution  matrix. 

This  defines  one  complete  cycle  of  the  use  of  ADAPS.  Each  part  can  be 
used  independently  of  the  others,  for  different  needs.  The  main  program  itself 
is  largely  at  the  discretion  of  the  user,  to  be  organized  as  best  suits  a  parti¬ 
cular  analysis  problem. 
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SECTION  III 


DEVELOPMENT  OF  A  NONLINEAR  DYNAMICAL  MODEL 
FOR  RIGID-BODY  MOTIONS 


The  ADAPS  is  a  six-degree-of-freedom  optimal  delivery  control  and 
aircraft-weapon  simulation  programming  system.  It  is  based  on  the  THRUST 
program  developed  in  Honeywell  [12]. 

The  equations  of  motion  in  ADAPS  are  referenced  to  a  flat  nonrotating 
earth.  The  forces  and  moments  acting  on  the  vehicle  and  the  weapons  a~e 
generated  by  gravity,  aerodynamic  effects,  and  thrust.  The  aerody 
coefficients  are  input  in  tabular  form  as  functions  of  Mach  number,  anb_e 
of  attack,  sideslip  angle,  etc.  The  aerodynamic  force  and  moment  sub¬ 
routines  allow  the  description  of  the  vehicle  and  weapon  characteristics 
over  a  wide  range  of  flight  conditions. 

ADAPS  has  three  running  phases:  (I)  modeling  and  linearization  phase; 
(II)  controller  optimization  phase;  and  (III)  performance  evaluation  phase. 
Briefly,  phase  I  consists  of  generating  a  prescribed  steady  trajectory  of  an 
aircraft  weapon  system  and  its  variational  equations  starting  from  a  given 
initial  condition  to  a  weapon  release  point,  then  continuing  to  generate  six- 
degree-of-freedom  free-fall  trajectories  of  a  weapon  and  its  variational 
equations  until  a  prescribed  target  altitude  is  reached.  The  coefficients 
of  the  variational  equations  are  obtained  by  a  simple  numerical  differentia¬ 
tion  process. 

In  the  following,  analyses  pertaining  to  phase  I  of  ADAPS  operation  are 
given. 


DEVELOPMENT  OF  THE  DIFFERENTIAL  EQUATIONS  OF  MOTION 

For  completeness,  this  subsection  presents  the  derivation  of  the  equation 
of  motion  of  an  airplane  and  a  weapon  [  13,  14,  15,  16,  17].  These  equations 
of  motion  are  implemented  in  subroutine  DYNK.  Since  they  are  unaffected  by 
the  interchangeable  subroutines  describing  alternate  airframe  and  weapon 
aerodynamics,  they  are  applicable  to  both  aircraft  and  weapon.  In  the  follow¬ 
ing,  aircraft  or  weapon  will  be  referred  to  as  a  body  or  a  rigid  body. 


Reference  Frames 

The  coordinates  necessary  to  specify  the  six  degrees  of  freedom  of  a 
rigid  body  are  defined  by  means  of  two  right-handed  reference  frames,  the 
x  ,  y  ,  z  frame,  or  earth-fixed  frame,  and  thex,  y,  z  frame,  or  moving 
frame  or  body  frame.  They  are  defined  as  follows: 
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•  Earth-Fixed  Frame  — 


0  =  origin,  fixed  at  the  target 

xe-axis  is  horizontal,  in  the  vertical  plane  containing  the 
initial  velocity  vector  of  the  mass  center  (downrange) 

ye-axis  is  crossrange  to  the  right 

ze-axis  is  vertically  downwards 

•  ■  Moving  Frame  — 

0  =  origin  at  the  center  of  gravity  of  body 

x-axis  is  parallel  to  the  longitudinal  axis  of  the  body 

y-axis  is  perpendicular  to  the  plane  of  symmetry  (for  a 
weapon  which  has  more  than  one  plane  of  symmetry 
choose  the  one  most  nearly  parallel  to  the  aircraft 
plane  of  symmetry  before  release) 

z-axis  is  perpendicular  to  xy  plane  and  positive  down 
viewed  by  the  pilot  (before  weapon  release  in  the  case 
of  the  weapon) 

Earth-fixed  and  moving  frame  definitions  for  aivplane  and  weapon  are 
illustrated  in  Figure  6. 


Figure  6.  Earth- Fixed  and  Moving  Coordinate  Frames 
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It  should  be  noted  that  the  orientation  of  earth-fixed  frame  differs  in  accor¬ 
dance  with  its  use  [171  The  earth-fixed  frame  generally  used  in  bombing 
and  range  tables  is  obtained  by  rotating  the  fixed  frame  described  here  by 
90  degrees  counterclockwise  about  the  x-axis. 

The  differential  equations  of  motion  consist  of: 

The  dynamical  laws  of  motion 

The  kinematical  relationships  between  the  reference  frames 


The  Dynamics  of  Motion 

Two  vector  equations  define  the  motion  of  a  rigid  body  completely.  They 
are  [13,  16,  17] relative  to  the  earth-fixed  frame. 

dvo 

F=ni— -  (3.1) 


where 

^  =  EfV  resultant  external  force  acting  upon  the  body 
vq=  velocity  of  mass  center  relative  to  inertial  frame 
m  =  ESmtotal  mass  of  the  body 

-4  —4 

T  =  ST.  resultant  external  torque  about  mass  center 
h  =  angular  momentum  of  the  body  about  mass  center 

For  a  rigid  body,  the  angular  momentum  about  mass  center  is  defined  as 


h  =  Sr  X  v  6m 

with 


(3.3) 


v  =  v  +  ui  x  r  (3.  4) 

o 

where 

v  =  velocity  of  the  mass  element  6m  relative  to  earth-fixed  frame 
v  =  velocity  of  center  of  mass  relative  to  earch-fixed  frame 
ou  =  angular  velocity  of  the  body  relative  to  earth-fixed  frame 
r  =  position  vector  from  c.g.  to  the  mass  element  6m 
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Substituting  <3. 4)  into  (3.3)  and  using  body- coordinates  one  obtains 


-♦-♦  2  2  2 

h  =  0)  £(x  +  y  +  z  )6m  -  £r(px+qy+rz)6m 


in  matrix  notation  further  simplification  yields  h  =  J  q 

y 


v’  -jre  J  is  the  moment  of  inertia  matrix  given  by 


I  -I  -I 

xx  xy  xz 

J  -  -I  I  -I 

xy  yy  yz 

\  *XK  ~VZ  *zz 


(3.5) 


In  order  to  supress  the  derivatives  of  moment  of  inertia  elements,  the 
equations  of  motion  given  by  (  3. 1)  and  (3.2)  are  expressed  in  body  axes 
which  rotate  with  the  body. 

Thus  (3.1  )  and  (3.2)  expressed  in  rotating  frame  of  reference  become 
F  =  m(g~j  +  m(iv  xvo)  (3.6) 


(3.6) 


T  =  |^+  5xh 
where  by  definition 

6t  dt  dt  dt 


(3.7) 


(3.7) 


Re  'rring  to  Equation  (3.7),  let 

ff  =  w  X  vq  =  T(qw-rv)  +  J(ru-pw>  +  S(pv-qu) 
In  matrix  notation  (3. 8)  becomes 


0  -r 
r  0 
-q  P 


')0 


(3.8) 


(3.9) 


Now  making  use  of  (3. 9)  equations  of  motion  (3. 6)  and  (3.  7)  can  be  expressed 


0  -r 


F  =m  V  +m 


(3. 10) 


P 

■J 


(3.  11) 


These  six  deferential  equations ,  called  the  Euler  equations  of  motion, 
completely  describe  the  dynamics  of  the  rigid  body. 


The  Kinematics  of  Motion 

To  describe  the  position  and  orientation  of  the  body  as  a  function  of 
time,  earth-fixed  frame  of  reference  is  introduced  as  defined  earlier.  To 
find  the  differential  equations  ot  the  coordinates  of  mass  center  relative  to 
the  fixed-frame  (i.  e.,  differential  equations  of  the  flight  path),  it  is  necessary 
to  develop  the  time  evolution  of  the  transformation  matrix  0(t)  which  trans¬ 
forms  a  vector  from  the  body  coordinates  to  the  earth-fixed  coordinates,  that 
is 

ve  *  fiv  (3. 12) 

Thus,  if  &e,  Jre,  fce  are  the  velocity  components  of  mass  center  in  the  earth- 
fixed  reference  system,  and  if  (u,  v,  w)  are  the  velocity  components  of  mass 
center  in  the  body  frame  it  follows  that 


=  0  v 


(3.13) 


Among  the  various  angular  coordinates  that  can  be  used  to  express  the 
direction  cosines  (i.  e.,  the  elements  of  fi  matrix),  a  system  of  Eulerian 
angles  is  commonly  used.  Because  of  their  trigonometrical  form,  how¬ 
ever,  Eulerian  angles  cause  a  singularity  in  the  differential  equations 
connecting  the  angular  velocity  with  the  Eulerian  angles  [173.  This  singu¬ 
larity  gi-;.^3  rise  to  a  considerable  tru?  nation  error  in  the  integration  pro¬ 
cess.  Tu  svoid  this,  the  components  of  a  normalized  quaternion  (versor, 
Euler  symmetrical  parameters)  ere  used  as  angular  coordinates.  For 
completeness,  discussions  on  both  the  usage  of  Eulerian  angles  and  quatern¬ 
ions  as  angular  coordinates  are  given  in  the  following  two  subsections,  re¬ 
spectively.  Quaternions  are  implemented  in  the  subroutine  DYNK  for 
nonlinear  simulation.  The  Eulerian  angles  are  used  for  the  perturbation 
model. 
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The  Eularian  Angles  as  Angular  Coordinates  --In  the  following,  the  elements 
of  the  transformation  matrix  Cl  are  expressed  in  terms  of  the  so-called 
Eulerian  angles.  In  general,  one  can  carry  out  the  transformation  from  a 
given  cartesian  coordinate  system  to  another  by  means  of  three  successive 
rotations  performed  in  a  specific  sequence.  The  Eulerian  angles  are  then 
defined  as  the  three  successive  angles  of  rotations.  Unfortunately,  there 
is  no  unanimity  in  the  literature  about  the  definition  of  the  Eulerian  angles. 

The  convention  of  C 17,  18],  and  [16]  will  be  adopted  here.  The  body  is 
imagined  first  to  be  oriented  so  that  its  axes  are  parallel  to  earth- fixed  axes. 
This  system  will  be  denoted  by  (Xj,  yj,  Zj).  The  sequence  will  be  started  by 
rotating  the  initial  system  of  axesjx,,  y^,  z.)  by  an  angle  if  clockwise  about 
the  z\  axis,  and  the  resultant  coordinate  system  will  be  labeled  the  xg,  yo*  *2 
axes.  In  the  second  stage  the  Intermediate  axes  Xg,  yg»  ^2  are  rotate^  about 
the  y*2  axis  clockwise  by  an  angle  B  to  produce  another  intermediate  set  Xg, 
yg,  zg.  Finally  the  xg,  yg,  Zg  axes  are  rotated  clockwise  by  an  angle  0  about 
the  xg  axis  to  produce  the  desired  x,  y,  z  body  system  of  axes.  Figure  7  illus¬ 
trates  the  various  stages  of  the  sequence.  The  Eulerian  angles  0,  if,  0  thus 
completely  specify  the  orientation  of  the  (x,  y,  z)  body  system  relative  to  the 
(x  ,  y  ,  z  )  earth-fixed  system. 

C  C  C 


Figure  7.  Body  Orientation  and  Euler  Angles 


The  elements  of  the  complete  transformation  E  can  be  obtained  by 
writing  the  matrix  as  the  triple  product  of  the  separate  rotations.  Let  Ej, 

Eg,  Eg  be  transformations  describing  the  rotations  t| i,  0,  and  0,  respectively. 
The  complete  transformation  E  from  earth-to-body  system  is  given  by 


where 

E  r  E,  E0  E 


Clearly,  rotation  about  Zj  axis  by  the  angle  *  is  given  by 
P cos  ijr  sin  ijf  0~1  / 


Ei  =  -sin  cos  (r  0 


L  0 


X1  “  xe 


Similarly  rotation  about  y2  axis  by  the  angle  6  is  given  by 


E2  = 


cos  6  0  -sin  0 


sin  0  0  cos  0 


Also  rotation  about  x3  axis  by  the  angle  0  is  given  by 


0  cos  0  sin  0 
0  -sin  0  cos  0 


The  product  matrix  E  then  follows  as 


cos0cosijr  cos0sin*  -sino 

E  =  -cos0sin$  +  sin0sin0cosifr  cos0cosi|f  +  sin0sin0sini|f  sin0cos0 

sin  0simji  +  cos0sin0cosi(r  -sin0cosi|r  +  cos0sin0sini|r  cos0cos0^ 

(3. 14) 

Since  E  is  an  orthogonal  matrix,  the  inverse  transformation  from  body  coor¬ 
dinates  to  earth- fixed  axes  is  then  given  by 


-sin0 


Q  =  E"1  =  E' 


(3.15) 
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Equations  (3. 13)  together  with  (3. 14)  and  (3. 15)  constitute  the  flight- 
path  equations  of  a  rigid  body  in  terms  of  body- axis  velocity  components 
which  in  turn  are  found  by  solving  Equations  (3. 10)  and  (3. 11).  We  develop 
in  the  following  the  differential  equations  of  the  Euler  angles  which  when  in¬ 
tegrated  enable  one  to  evaluate  the  direction  cosine  matrix  ft  as  a  function 
of  time. 

Let  T)  be  an  arbitrary  fixed  vector.  Let  r\e  and  t)  be  its  component 
representation  in  the  earth-fixed  and  body  frames  respectively. 


Then  one  can  write 


r\  =  E  T» 


(3.16) 


Differentiating  (3. 16)  w.  r.  t  time  and  noting  that 


yields 


tl  =  ET»e 

using  (3. 16)  in  (3. 17)  yields 
r\  -  E  Er  ri 


(3. 17) 


(3.18) 


where 


%  dTk 


'I  “  col|df  ’  “dt  ’  dt  j 

On  the  other  hand,  the  time  derivative  “n  w.  r.  t.  rotating  frame  of  reference 
is  given  by 


=  M.  +  Z  x  n 

dt  6t  +  u  X  T) 


where 


(3.19) 


*1*  r  2*  +-r  d_l  +s  dA 

1  +  3  Ht  +  K  rtf 


(3. 19a) 


In  matrix  notation  (3. 19)  and  3. 19a)  become 


T)  =  -W  T) 


(3.20) 
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When  q©  =  0,  the  quaternion  is  said  to  be  a  pure  quaterion.  By  definition  then,  a 
vector  is  a  pure  quaternion  and,  consequently,  the  bracketed  term  in  (3.25) 
is  called  the  vector  part  of  q.  Defining  four  angles  a,  3,  y,  and  cp  in  the 
following  way,  one  obtains  polar  representation  of  a  quaternion  [17] 

q  =  (qo  +  iqx  +  jqg  +  kq3)  =  VN(q)  [cos  +  sin  j  (i  cos  a  +  j  cos  p  +  k  cos  y)J 

(3.27) 


The  conjugate  of  q  is  defined  as 

q  =  %  ~  icli  ~  W2  '  k%  (3.28) 

It  may  be  verified  that 

qq  =  qq  =  N(q)  (3.29) 

An  equivalent  definition  expresses  the  quaternion  as  a  4  x  4  matrix  of  a  spe¬ 
cial  type 


q  «*  Q 


%  ql  q2  q3 

-ql  %  "q3  q2 

-q2  %  %  -qx 

-q3  -q2  ql  qo 


(3.30) 


where  q_,  q^,  q2,  q^  are  the  components  of  a  quaternion  given  in  (3.25)  .  In 
matrix  rorm  the  conjugate,  q,  of  the  quaternion  q  becomes  Q the  transpose 
of  Q. 

q  -  Q  '  (3.31) 

The  j  roduct  r  =  qp  of  two  quaternions  represents  an  operation  on  p  that  rotates 
p  by  an  amount  cp/2  and  stretches  its  magnitude  by  VN(q).  Let 

r  =  rQ  +  irx  +  jr2  +  kr3 

q  =  %  ♦frqj  +  jqx  +  kq3]  (3.32) 

P  =  P0  +  iPj  +  jp2  +  kp3 

Using  (3.26  )  and  (3. 26a)  the  product  r  *  qp  can  be  expressed  as 
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(3.33) 


ro  "  (qopo  '  qlpl  ‘  q2p2  •  q3p3> 

rl  =  (qlpo  +  qopl  •  q3p2  +  q2p3> 

r2  =  <q2po  +  q3pl  +  qop2  ‘  qlp3> 

r3  =  <Vo  •q2pl  +  qlp2+qop3) 

The  same  result  can  also  be  obtained  by  considering  the  matrix  equivalent  of 
r  =  qp, 

R  =  QP  (3.34) 

and  equating  elements  of  the  first  row  to  find  r  ,  r ^ ,  r^,  r^ 

Clearly  (3.33)  can  be  written  as 

r  =  P'q  (3.35) 

where  r  and  q  now  are  4x1  column  vectors.  Equation  (3.  33)  can  also  be 
written  as 

r  =  Qp  (3.36) 

where  p  is  the  4x  1  vector  repres entation of  the. quaternion  and Q  is  a  4x4  ma¬ 
trix  obtained  from  Q  by  interchanging  the  first  row  and  the  first  column  of  Q- 

Besides  rotating  and  stretching  four- dimensional  vectors,  quaternions 
also  effect  rotations  and  stretchings  of  the  three-dimensional  space,  and  it 
is  this  property  that  makes  quaternions  useful  for  the  description  of  rigid 
rotations.  To  demonstrate  this,  let  x  be  a  pure  quaternion  (i.  e. ,  a  vector) 
described  by 

x  =  +  jx2  +  kx3  (3.37) 

Let  X  be  a  unit  quaternion  (that  is,  XX  =1)  described  by 

X  =  \5  +  iX1  +  jx2  +kX3  (3.38) 

and  Xbe  the  conjugate  of  X. 

Now  consider  the  quaternion  y  defined  by  the  triple  product 
y  =  Xx  *  (3.39) 
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By  actually  performing  the  indicated  multiplications  with  the  rules  given  in 
(3.26)  and  (3.26a),  or  writing  (3.39)  with  the  aid  of  (3.35)  and  (3.36)  as 


r\  =  kc  ♦*  Ax 
y  =  (\x)  X  =  rjX 


(3.39a) 

(3.39b) 


and  substituting  (3,39a)  into  (3.39b)  upon  returning  from  four  to  three 
dimensions  for  y  one  obtains  the  equation 

y  =  flx  (3.40) 


where 


Q 


2U02+\12)-i  attjWa)-  2(xix3+xox2> 
2(\j\2+^0X3^  2(Xo  +x2  ^  *  2^X2x3"XoXl^ 


ZfXjXg-X^g)  2(X2X3+XoXl* 


2<Xo2+X32)-  *> 


(3.41) 


It  can  readily  be  shown  that  (3.39)  preserves  the  norm,  that  is 
xx  =  yy 


(3.42) 


On  the  basis  of  (3.40)  and  (3.42)  one  concludes  that  (3.  39)  defines  a  rotation 
of  x.  This  demonstrates  that  quaternions  can  be  used  as  angular  coordinates. 
Now,  the  differential  equations  of  the  elements  of  a  quaternion  will  be  developed. 

Let  ?  and  t)  be  fixed-vectors  attached  to  earth-fixed  frame  and  rotating 
body  frame  respectively  as  shown  in  Figure  8. 


Figure  8.  Relations  between  Fixed  and  Rotating  Frames 
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Consider  transformation  defined  by  (see  Equation  3.  39) 

r\  =  MX  (3.43) 

To  keep  track  of  reference  frames,  rewrite  Equation  (3. 43)  in  the  body 
frame  as: 

=  X  ?bX  (3. 43a) 

where  the  subscript  b  implies  the  corresponding  vectors  expressed  in  body 
frame  (x,y,x).  Let  us  further  assume  that  the  two  vector  systems  coincide 
when  X  =  1. 

That  is 

\  =  ?e  (3‘44) 

where  the  subscript  e  implies  the  corresponding  vector  expressed  in  the 
earth-fixed  frame.  Since  r\  moves  with  the  body  frame;  (3.44)  holds  for  all 
values  of  X. 

Substituting  (3. 44)  into  (3.43a)  yields 

?e  'X^X  (3.45) 

or 

?b*X?eX  (3.45a) 

-4 

This  relates  in  quaternion  notation  the  components  of  the  vector  5  in  the 
two  coordinate  frames.  In  matrix  notation  it  becomes 

Y  =  A'  Y  A  (3.46) 

e 

where  Y  is  the  matrix  representation  of  the  quaternion  ?  . 

Differentiating  this  and  noting  that  Yg  is  a  fixed  quaternion  matrix  yields 

Y  =  *A#Y  A+A'Y  A  (3.47) 

e  e 

but  from  (3.  46) 

Ye  =  AYA'  (3.48) 
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Substituting  this  into  (3.47)  and  simplifying  yields 

Y  =  A'AY  +  YA'A 

Noting  that  y  is  a  pure  quaternion  so  that  [see  Equation  (3.  30)1  Y'=  -Y 
results  in 

Y  =  A'  A  Y  •  (A/AY)/  (3.49) 

Now,  differentiation  of  a  vector  in  a  rotating  frame  in  vector  notation  gives 
Equation  (3. 19): 

?  =  UX?  (3,50) 

which  can  be  shown  to  have  the  quaternion  representation 

Y  =  (“  2  )  CWY  -(WY )#]  (3.51) 

where  W  is  the  quaternion  matrix  of  the  angular  velocity  vector  of  the  rota¬ 
ting  frame  presented  in  the  rotating  *>ame,  that  is 

5J  =  Ip  +  jq  +  Kr  (3.  52) 

Comparing  (3.49)  and  (3.51)  gives 

A 'A  =  W 

or 

A'=<-|0WA'  (3i63) 

In  the  manner  of  Equation  (3.  34),  after  performing  the  multiplication  the  first 
column  of  this  equation  constitutes  the  differential  equations  of  the  components 
of  the  quaternion  which  are  given  by  [see  Equation  (3.30)3 


To  integrate  this  set  of  equations,  an  initial  value  of  X  must  be  specified. 


Initialization  of  the  Angular  Coordinates  —  By  comparing  corresponding 
elements  of  matrices  (3. 14)  and  (3.41),  the  components  of  the  vector 
can  be  related  to  Euler  angles  [21]  as  follows 


Xo 

ill 

=  cos 

cos 

1  cos 

<t>  . 
I  + 

sin  | 

sin  |  sin  4 

=  cos| 

cos 

-fsin 

1  . 
2 

sin| 

sin  j-  cos  4 

X2 

=  cos  | 

sin 

f  cos 

!♦ 

sin  | 

cos  4  sin  4 

X3 

=  “  cos  l 

sin 

e  . 

?  sin 

l  + 

sin  | 

6  0 
COS  Y  cos  \ 

(3.55) 


This  set  of  equations  can  be  used  to  obtain  initial  values  of  the  versor  (unit 
vector)  components.  An  alternate  way  of  obtaining  the  versor  components 
from  the  Euler  angles  is  to  compute  the  initial  cosine  matrix  and  use  the 
following  easily  verified  relations: 


2  1-ftrE 
Xo  '  4 

where  tr  is  the  trace  operator,  and 


provided  that  XD  ^  0.  Since  for  certain  initial  conditions  (i.  e. ,  8  -  tt,  \b  =  0 
<p  -  u)  this  condition  is  violated.  Equation  (3.  55)  is  complemented  in  ADAPS 
instead  of  (3.  56) 


Gravity  Components  in  the  Body  Axes 

The  inertial  components  of  gravity  are  resolved  into  the  required  body 
components  by  using  the  direction  cosine  matrix. 


which  picks  the  third  column  of  E  given  in  (  3. 14),  that  is 


35 


■ 


:\»\7  ****-:.  t 


r. ,*=■•-* tvs?^  V*v> »  .*  .  ; 


L  \ 


8x1 

ey 

=  g 

V  *1 

\ 

-sin  8 
cos  8  sin  0 
cos  0  cos  0 


\ 

/ 


This  finishes  the  analysis  of  the  equations  of  motion. 


(3.59) 


Summary  of  the  Analysis 

In  the  development  of  the  equations  of  motion,  the  aerodynamic  forces 
are  taken  to  be  through  the  aircraft  mass  center  and  the  moments  are  about 
the  body  axes.  The  accelerations  on  the  aircraft  are  computed  by  combining 
the  aerodynamic  forces  and  the  forces  from  the  engine.  All  cross  products 
of  inertia  are  included  to  allow  for  a  nonsymmetricel  body. 


Accelerations  due  to  Aero  and  Thrust  Forces  in  Body  Coordinates 


V 

_  -L 

m 

”x  +  i:~ 

a 

y 

a„ 

Y 

z  +  zm 

z 

L  TJ 

(3.60) 


Total  Moments  in  Body  Coordinates 


M 

L  +  Lrp 

X 

T 

My 

= 

m+mt 

Mz_ 

N+  Nt 

(3.61) 


*  Differential  Equations  of  Body  Translation  Velocities  in  Body 
Coordinates  — 
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with  prescribed  initial  conditions 


•  Differential  Equations  of  Body  Angular  Velocities  in  Body 
Coordinates  — 


is  a  symmetric  matrix;  or 


where  the  coefficient  matrix  D  is  J_1and  symmetric  and 


(3.63) 


(3.64) 


(3.65) 
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and 


■TT^B^Tr*5^: 


I 


p 


K  : 


The  Differential  Equations  of  Body  Mass  Center  in  the  Earth-Fixed 
Coordinates  U.  e. ,  Flight  Path)  - 


=  E' 


u 

v 

w 


(3.72) 


with  prescribed  initial  conditions 


*e«» 


ye<o) 

ze(0) 


DEVELOPMENT  OF  INTEGRATION  ALGORITHM  FOR  THE  EQUATIONS  OF 
MOTION 

The  method  of  integrating  the  differential  equations  of  motion  in  ADAPS  is 
the  open  quadrature  process  (i.  e. ,  Adams  formula)  0j2].  It  is  given  by 

xk  =  xk-i  +  h(1+T  Vi  (3.73) 


where 

h  =  He  ‘  tk-: 


x  =  f(x,  t) 


(3.74) 

(3.75) 


and 


v  =  backward  difference  operator 

Truncation  error  associated  with  the  process  is  of  the  order  of  h3  and  is  given 

by 

E  =  a2h3x(3)(S),  tk_2  <  tfc  (3.76) 

This  integration  process  can  be  written  as 

H - t - h- 

V2  Vl  *k  <3‘7?) 


xk  =  Vl+l(3Vl  "  *xk-2) 


To  start  the  solution,  the  derivatives  at  two  time  points  are  needed.  The 
solution  can  be  started  by  using 

Y  S  V  X  Vi  vr 


tk-l<5< 


(3.79) 


with 


E  = 


h2  x  (§) 


X 


In  ADAPS,  the  initial  values,  x.j,  are  set  to  zero  for  simplicity.  For 
small  h  the  error  caused  by  this  simplification  is  small. 
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SECTION  rv 


DEVELOPMENT  OF  TOTAL  FORCE  AND  MOMENT  SYSTEM 
FOR  RIGID  BODIES  MOVING  IN  UNSTEADY  AIRMASS 


The  methods  by  which  the  aerodynamic  and  thrust  forces  and  moments 
are  introduced  into  the  six-degree-of -freedom  trajectory  program  are  pre¬ 
sented  in  this  section.  Aerodynamic  forces  and  moments  for  aircraft  are 
treated  first.  The  computer  program  which  implements  the  aerodynamical 
model  of  aircraft  is  called  subroutine  AERK.  The  aerodynamic  -  of  bombs 
are  treated  next.  The  corresponding  subroutine  is  called  subroutine  WAERK. 
Then  a  model  for  the  thrust  forces  and  moments  is  developed.  The  program 
which  implements  this  model  is  called  subroutine  THRUSK.  Finally,  the 
effects  of  moving  air  mass  on  the  aerodynamics  of  rigid  bodies  are  treated 
to  take  into  account  the  mean  winds  and  stochastic  wind  gusts.  The  program 
which  implements  the  wind  model  is  called  subroutine  WINDK. 


AERODYNAMIC  FORCES  AND  MOMENTS  FOR  AIRCRAFT 

The  aerodynamic  fo  ces  and  moments  are  computed  by  making  use  of 
extensive  aerodynamic  data  tables.  The  primary  objective  of  the  function 
lookup  subroutine  (FLOOR),  presented  in  Volume  II  is  to  provide  for  a 
complete  accounting  of  the  various  contributions  to  the  aerodynamic  forces 
and  moments  regardless  of  the  flight  conditions  or  the  body  (i.  e. ,  aircraft, 
weapon)  being  considered.  The  technique  used  is  an  n-dimensional  table 
lookup  and  linear  interpolation.  This  method  has  the  advantage  of  accurately 
describing  even  the  most  nonlinear  variations  with  a  minimum  of  preparation 
effort.  However,  the  amount  of  storage  space  and  computing  time  increases 
rapidly  with  the  number  of  dimensions  of  the  tables. 


The  Functional  Form  of  the  Aerodynamic  Forces  and  Moments 

The  major  variables  that  affect  the  aerodynamic  characteristics  of  a 
body  are:  Mach  number,  M;  dynamic  pressure,  q;  angle  of  attack,  a;  angle 
of  sideslip,  8;  total  linear  velocity,  angular  velocity,  and  control  surface 
deflections.  Two  typical  functional-dependence  relations  can  be  written  as 

F  =  F(Xj,  Xj,  Xg,  x2»  h*  M)  (4.  1) 

=  F(V,  a,  8,  V,  a,  B,  Xg,  Xg,  6,  6,  h,  M)  (4.  2) 
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where 

Xj  =  linear  velocity  state  vector 

x2  =  angular  velocity  state  vector 

6  =  control  surface  deflection  vector 

h  =  altitude 

M  =  Mach  number 

V  =  total  velocity 

a  =  angle  of  attack 

p  =  sideslip  angle 

The  decomposition  of  (4. 1)  or  (4.  2)  into  functional  relationships  with 
fewer  number  of  arguments  is  needed  for  practical  reasons.  Unfortunately, 
there  is  no  standard  form  for  this  decomposition,  and  it  is  dependent  on  the 
available  data.  Usually  airframe  and  weapon  manufacturers  provide  empiri¬ 
cal  or  estimated  relations  in  the  form  of  curves  or  data  tables. 


Usage  of  Aerodynamic  Data 

Aerodynamic  characteristics  strongly  depend  upon  the  orientation  of  the 
relative  wind  or  velocity  vector  with  respect  to  body,  so  the  angle  of  attack, 
a,  and  sideslip  angle,  0,  which  define  the  orientation  with  respect  to  the  air 
mass  are  used  often,  as  indicated  in  equation  (4. 2).  In  general,  the  aero¬ 
dynamic  data  is  classified  into  two  groups;  static  data  and  dynamic  data. 
Static  data  implies  that,  during  the  wind  tunnel  testing,  the  body  is  at  rest 
with  respect  to  the  relative  wind.  The  dynamic  data  implies  that  body  oscil¬ 
lates  or  rotates  with  respect  to  the  relative  wind. 

The  aerodynamic  data  are  usually  given  in  the  "wind-tunnel  stability 
axes"  defined  as  follows  C 15] ; 

Origin:  center  of  mass  (eg) 

x  :  in  the  direction  of  motion  along  the  projection  of  V  upon 

the  body  reference  plane. 

y_:  same  as  body  axes,  positive  right 

s 

z  :  perpendicular  to  x  y  plane  forming  right-hand  triad 

s  s  s 


The  wind  tunnel  stability  axes  system  is  illustrated  in  Figure  9.  Clearly 
the  body  axes  system  and  stability  axes  systems  are  related  to  each  other  by 


where 


The  resolution  of  the  total  aerodynamic  force  in  the  xz  plane  is  shown  in 
Figure  10  where  lift,  L,  and  drag^D,  are  forces  normal  and  parallel, 
respectively,  to  the  projection  of  V  in  the  xz  plane.  Lift  and  drag  are  defined 
to  be  positive  as  illustrated. 


Aerodynamic  forces  and  moments  are  expressed  in  terms  of  the  basic 
aerodynamic  coefficients  in  the  wind  tunnel  stability  axes  with  origin  located 
at  an  arbitrary  reference  point.  The  aerodynamic  force  and  moment 
coefficients  are  defined  by  the  following  relations: 


(4.4) 


where 

_  o 

q  =  dynamic  pressure  (lb/ft  ) 

2 

S  «  wing  area  (ft  ) 

b  =  wing  span  (ft) 

c  =  wing  mean  aerodynamic  chord  (ft) 


44 


It  is  assumed  that  the  moment  reference  center  0ca  is  located  by  Arca 
from  the  eg  of  the  aircraft  as  shown  in  Figure  11. 


Figure  11.  Moment  Reference  Center  With 
Respect  to  Mass  Center 


The  total  aerodynamic  forces  at  the  eg  in  body  axes  are  then  given  by 


f  =  E'  f  **  E'  Y 
a  s  sa  s  s 


(4.6) 


Similarly  the  total  aerodynamic  moments  at  the  eg  about  the  wind  tunnel 
stability  axes  are  given  by 


=  M  +<Ar  )  x  f 


sa  sea 


(4.7) 


Or  in  matrix  notation,  and  in  body  axes 


M  =  E'  Mo  +  AR  f 
a  s  sea  ca  a 


(4.8) 


where,  in  terms  of  a  reference  point  RP  on  the  body  (see  Figure  11), 


Ar  =  r  +  r 
ca  eg  ca 


(4.9) 


Ar  „  B  position  vector  from  0  to  0 
ca  r  eg  ca 

»  position  vector  from  0  „  to  RP 
eg  r  eg 

r„„  *  position  vector  from  RP  to  0 _ 

ca  r  ca 


and,  with  coordinates  expressed  in  body  axes. 


Substituting  (4.  5)  and  (4. 9)  into  (4.  8)  yields  the  moment  components  in 
body  axes  at  the  eg. 


(4.11) 


Total  Aerodynamic  Coefficient  Model  for  Aircraft 

In  the  formulae  given  below  the  superscript  "oH  indicates  degrees,  the 
subscript  "a"  denotes  quantities  with  respect  to  the  air-mass,  and  subscript 
"w"  stands  for  wings  with  respect  to  the  air  mass  (aw  =  aa  +  iw).  This 
notion  will  be  useful  in  treating  a  moving  air  mass  as  discussed  later.  The 
aerodynamic  force  coefficients  in  the  wind  tunnel  stability  axes  are  assumed 
to  be  in  the  following  form  [36]; 
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C  (h,M  )  C  <h-VMa> 
6a  I  6s 


+  cy6a<Ma>  0  Cy5r<h-Ma>|  6”s 
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*«a  a 
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CDS  <«°w>IF0bp 
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eLsp(h,Ma)  vr  (a°wTira)-c—  (a«p  6°lg 


Similarly  the  aerodynamic  moment  coefficients  in  the  stability  axes  are 
assumed  to  be  in  the  form  of 


Cmsc<Ma-h'a  w> 


0  CjfaVMy  o  0  0 

+  0  0  I  ao  +  Cm.<h-Ma>  0 
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Cn<h-“°WMa>  0  Cn  (Ma'h-a°w)  0  0  IST  ra 
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6a  6r 


AERODYNAMIC  FORCE  AND  MOMENT  MODEL  FOR  BOMB 

A  complete  aerodynamic  model  for  a  slowly  spinning,  four-finned  bomb 
is  given  in  [17],  where  the  aerodynamic  parameters  are  assumed  to  be 
linear  functions  of  spin,  cross-spin  and  accidental  configurational  asymmetry 
but  nonlinear  functions  of  yaw  orientation  and  roll  orientation.  Then  the 
effects  of  roll  orientation  on  the  aerodynamic  forces  and  moments  are 
obtained  by  a  Fourier  series  expansion  of  roll  angle  (the  angle  between  the 
[plane  of  yaw]  and  a  reference  fin).  The  model  based  on  [  17]  (i.  e. ,  Cohen's 
model)  requires  approximately  20  aerodynamic  tables  (i.  e. ,  tests),  and  these 
tables  are  not  readily  available. 

In  ADAPS,  a  simplified  aerodynamic  bomb  model  is  developed.  It 
utilizes  generally  available  bomb  aerodynamic  data.  The  effect  of  roll 
orientation  is  ignored.  [The  cross -velocity  frame  and  cross-spin  frame  are 
the  same  as  defined  in  [17].] 


Simplified  Aerodynamic  Model  for  Bomb 

The  reference  axes  used  in  the  bomb  aerodynamic  model  are  illustrated 
in  Figure  12. 

0  *  origin,  at'  the  [  center  of  gravity]  of  bomb 

The  weapon  body  axes  are  defined  as  follows: 

x-axis  is  along  the  bomb  body,  positive  forward 

y-axis  is  horizontal  positive  right 

z-axis  is  perpendicular  to  the  xy  plane,  positive  down. 

The  cross-velocity  axes  are  defined  as  follows: 

x.  is  the  same  as  the  body  x-axis 

*  -*222 
z.  is  in  the  direction  of  cross- velocity  v  (V  *  v  +  w  ) 

X  Cel  Ca 

yj  is  perpendicular  to  the  XjZ^  plane  forming  a  right-handed  system. 
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Figui  e  12.  Definition  of  Stability- Like  Axes  for  Bomb 


The  cross-spin  axes  are  defined  as  follows: 

x„  is  the  same  as  the  body- axis 

•>  _  2  2  2 

y2  is  in  the  direction  of  cross-spin,  qc<qc  =  q  +  r  ) 

z2  is  perpendicular  to  x2y2,  the  plane  forming  a  right-handed  system. 
The  aerodynamic  forces  and  moments  in  body  axes  are  given  as 


=  qS  SC 


M  =  qSd  EC 


I  " 
!  1 


where 


q  =  1/2  p  vf.dynamic  pressure,  lb/ft2 

a 

d2  2 

S  =  cross-sectional  area  of  the' bomb  =  tt  -j— ,  ft 
d  =  diameter  of  the  bomb,  ft 

and  (EC,  EC  ,  EC)  -  aerodynamic  force  coefficients  in  body  axes  system 
x  y  z 

(ECj,  ECm,  SCn)  =  aerodynamic  moment  coefficients  in  body  axes  system 

The  nondimens ional  aerodynamic  data  are  given  in  the  cross-velocity 
frame.  They  are: 

C^(M)  =  axial  force  coefficient,  along  x-axis,  positive  aft 

Cjj(a,  M)  =  normal  force  coefficient  perpendicular  to  the  plane  of  yaw 

a 

CN^(a,  M)  =  coefficient  of  normal  force  due  to  canted  fin  shielding 
* 

Cm(a#M)  =  restoring  moment  coefficient,  about  y^^ -axis 

Cms (a,  M)  =  restoring  moment  coefficient,  due  to  canted  fin  shielding 
about  y ^  -  axis 

Cmq(a,  M)  =  damping  moment  coefficient  about  y2~ axis 

The  data  for  the  moment  coefficients  are  referenced  to  the  center  of 
gravity  [26,  29l 

The  transformation  matrices  from  the  cross -velocity  frame  to  body  frame, 
and  from  the  cross-spin  frame  to  body  frame  are  given  respectively  as 

*  9  Ei*i  n  =  E2  T| 


where 


and 

V 

0,  =  tan"1  —  and  0,  =  tan"1  —  ,  deg 
1  wa  2  q 

The  aerodynamic  coefficients  along  the  body  axes  are  given  by 


0  cos  0^  sin 
0  -sin  cos  <t>^ 


0  cos  02  sin  Pg 
0  -sin  <f>2  cos  02 


m 


11 


sc 


( 


-Ca(M) 


sc 


=  E, 


sc. 


-CN(a,  M)  -C*T(a,  M)  6 


N 


6 


/  lci\ 


l 


SC. 


m 


\N 

where 


C  (a,  M)  +  C  (a,  m)6  +E0  [C  (a,.M)  2V 


m 


m 


m. 


V 


a  =  tan 


-1  vca 


,  magnitude  of  yaw  (deg) 


ca 


a 


=  ~\jv^  +  w a2  ,  cross -velocity  (ft/sec) 

Vu  2  +  v  2  +  w  2  ,  total  velocity  (ft/sec) 

a  a  a 


V2  2 

6z  +  6  y  ’  ma6nitude  *in  can^  (de6) 


V2  2 

q  +  r  ,  cross-spin  (deg/sec) 


DEVELOPMENT  OF  A  MODEL  FOR  THRUST  FORCES  AND  MOMENTS 


In  general,  the  total  thrust  forces  and  moments  acting  on  a  rigid  body 
depend  upon  the  positions  and  orientations  of  the  thrust  producers  with  respect 
to  the  body  axes  ana  the  magnitudes  of  the  thrusts  (geometry).  There  are 
also  dynamics  associated  with  the  thrust  variables  since  they  are  produced 
and  oriented  by  engines  and  actuators. 


To  provide  an  analysis  tool  by  which  the  effects  of  various  control  points 
and  methods  can  be  investigated,  both  aspects  are  considered  in  the  develop¬ 
ment  of  a  thrust  model  in  ADAPS. 


In  the  following,  a  geometric  model  for  the  effective  thrust  acting  on  a 
rigid  body  is  developed  first.  Then  the  dynamics  of  thrust  producers  are 
treated.  In  the  development,  the  effects  of  angular  momentum  of  rotating 
thrust  producers  are  neglected. 
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In  the  following,  first  the  description  of  a  geometrical  model  for  thrust 
producers  is  given.  Then  the  total  force  and  moment  contributions  of  thrust 
producers  are  developed  in  the  form  of 


fT  =  Bt  (xd)yT  (4. 12) 

mT  =  Bj,  (xd)yT  (4. 13) 

where 


fT  =  total  thrust  force  vector  along  body  axes 
nirp  =  total  thrust  moment  vector  along  body  aves 
Brj,  =  thrust  force  coefficient  matrix  of  size  3  x  v 

A 

B,p  =  thrust  moment  coefficient  matrix  of  size  3  x  v 

x,  =  state  vector  of  thrust  orientation  actuator  positions 
d 

y>P  =  effective  thrust  output  magnitude  vector 
v  =  number  of  thrust  points  in  the  thrust-producing  system 
Figure  13  shows  the  geometry  of  a  single  thrust  producer. 
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and 


r*  Tj 
SrTr  4yTj 

Ki 


«  g  J 


(4.  20) 


Let  Tj  be  the  effective  thrust  magnitude  of  the  jth  thruster.  Then,  in  terms  of 
of  the  elevation  and  azimuth  angles,  ©tj  and  \|fTj  of  the  nozzle  axis  from  the 
body  axis,  the  components  of  thrust  forces  fppj  along  body  axes  can  be  ex¬ 
pressed  as  J 


f„.  =  b.T. 

Tl  3  3 


3  =  1,  2,  ..  ,v 


(4.21) 


where 


cos  *Tj  COS  6Tj  I 
sin  ^ 

-cos  ^  sin  @T_. 


Making  use  of  \4.  21 )  and  (4. 22)  yields 

f«P  =  Brj.y.j, 


(4. 22) 


(4.  23) 


where 


bt  =  [bilb2l"-lbJandy 


=  T 
T  *2 


Similarly,  using  (4. 21)  in  (4, 18)  yields 


mT  =  Bt  uT 


(4.  24) 


(4.25) 


where 


BT  ^  bl  I  b2  I  **•  I 


(4.  26) 


£ 


with 


b.  =  ARm.b.  and 


(4.  27) 


"(AyTj)(’c3j)  -  (AZrpjKbgj) 

(AzTjWb1j)  ”  <AxTj)(b3j) 

(^x^xbgj)  -  (^yTj)(b1j) 


(4.  28) 


Equations  (4.23)  and  (4.25)  constitute  the  geometry  of  the  thrust  producers. 
In  the  following,  the  dynamics  of  thrust  producers  and  thrust  orientation 
actuators  are  treated  briefly. 

Dynamical  Model  for  Thrust  Magnitudes 


It  is  assumed  that  the  magnitude  dynamics  of  each  thrust  producer  can  be 
represented  by  a  first-order  transfer  function.  This  is  shown  in  Figure  14. 


Figure  14.  Dynamical  Model  for  a 
Thrust  Producer 


In  this  model  <M 

fVi 

x,p(j)  =  the  magnitude  state  of  the  j  n  thrust  producer  in  percent 

yT(j)  =  effective  thrust  output  of  the  j^1  thrust  producer  acting  on  the 

rigid  body  in  lbs 

h^,(j)  =  nonlinear  output  function 

til 

Ufp(j)  =  throttle  input  to  j  thrust  producer  in  percent 
a  (j),  b^fj)  =  transition  and  input  coefficients  of  the  producer 


^•"iyi-H 


The  output  function  for  the  two  main  thrust  producers  (i.  e. ,  engines)  is  assumed 
to  be  in  the  following  form  [36]: 

yT(j)  =  C(l+e1M)(l+e0h)]yo  +  |(l+e2M)(l+e0h)(c)[xT(j)-§T(j)]| ,  j  =  1,  2 

In  this  equation  eg,  e,  and  e  2  are  Mach-number-  and  altitude- dependent  co¬ 
efficients  of  the  effective  thrust,  yQ  is  thrust  bias,  c  is  the  thrust  output 
coefficient.and  ?ij<  (3)  is  a  function  of  Xpij).  These  coefficients  are  piece  wise- 
constant  functions  of  xT(j),  and  typical  values  for  F-4  engines  are  given  in 
Table  I. 

Dynamical  Model  for  Thrust  Orientation  Actuators 

The  thrust  orientation  state  of  each  thrust  producer  is  described  by  the 
azimuth  and  elevation  angles  of  the  nozzle  axis  as  defined  in  Figure  13. 


♦  (j) 

xd(j)  w 

Le  <j>] 


(4.  29) 


Table  I .  Effective  Thrust  Output  Coefficients,  j  =  1, 2 


[Coefficient 


0  *xT  *  50 


50  £  Xrj,  £  100 
9650 


-0.25150 

-<0.25)10~4 


-0.2342 
0.  32846 
-(0.  25)10-4 


It  is  assumed  that  first-order  dynamics  is  associated  with  thrust  orientation 
actuators.  Thus  for  each  thrust  producer  (see  Figure  14b) 

£  .  =  A  .V .  +  B  .u  .  (4,  30) 


=  A^x^  +  B^Uj 
^d  -  ^do  +  Xd 


(4.  31) 


where 


^do 


orientation  bias  * 


(4.  32) 


Implementation  of  the  thrust  orientation  dynamics  is  outside  the  scope 
of  the  present  program.  However,  in  ADAPS,  the  thrust  model  includes 
equation  (4.  31)  through  which  thrust  orientation  dynamics  can  be  inserted  into 
the  overall  program. 


EFFECT  OF  MOVING  AIR  MASS  ON  THE  AERODYNAMICS  OF  A  RIGID  BODY 

The  aerodynamic  force  and  moment  system  developed  above  is  with  respect 
to  the  air  mass  (i.  e. ,  atmosphere).  It  is  known  that  the  air  mass  through 
which  a  rigid  body  flies  or  Mis  is  in  a  motion  which  is  variable  both  in  time 
and  in  space  (Figure  15).  In  this  subsection  the  influence  of  this  motion  on 
the  rigid-body  aerodynamics  is  presented. 


Modeling  for  the  Influence  of  Unsteady  Air  M ass 

Meteorological  observations  indicate  that  the  velocity  field  of  air  mass 
(i.  e. ,  winds)  in  the  lower  atmosphere  consists  of  two  distinct  components, 
a  low-frequency  component  with  energies  concentrated  in  0. 01- cycle/hour 
range  and  a  high-frequency  component  with  energies  in  the  70-cycle/hour 
range.  The  former  is  called  the  "mean  wind",  and  the  latter  is  referred  to 
as  atmospheric  turbulence  or  simply  as  "wind  gust". 


Figure  15.  Rigid  Body  in  an  Unsteady  Air  Mass 
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In  the  earth  axes,  the  velocity  field  of  the  air  mass  can  thei«=*ore  be 
decomposed  into 

w  =  w  +  w  (4. 33) 

e  e  e 

where,  in  matrix  notation, 

we  =  velocity  field  vector 

wg  =  mean  velocity  field  vector 

we  =  gust  velocity  field  vector 

and  all  are,  in  general,  position-  and  time-dependent. 

The  influence  of  motion  of  the  air  mass  on  rigid-body  aerodynamics  is 
taken  into  account  with  various  degrees  of  accuracy  by  considering  a  rigid 
body  as  a  point,  line,  surface  or  volume.  When  a  rigid  body  is  considered 
as  a  point  in  the  air  mass,  then  the  relative  motion  appears  as  the  difference 

v  =  v  -  w  (4. 34) 

cl 

where  w  is  the  velocity  field  vector  in  body  axes.  In  this  case,  the  aero¬ 
dynamic  effect  of  the  motion  of  air  mass  is  accounted  for  by  the  use  of  the 
equivalent  (i.  e. ,  producing  the  same  aerodynamic  effects)  velocity  Va,  angle 
of  attack  a  and  angle  of  sideslip  0  . 

cL  a 


It  should  be  noted  that  the  equivalent  linear  velocities  are  to  be  used  only 
for  developing  the  aerodynamic  forces  and  moments.  Elsewhere  in  the 
dynamical  equations  the  inertial  velocities  are  used. 

If  a  rigid  body  is  assumed  to  have  one  or  more  dimensions  in  space,  then 
the  space  distribution  of  the  velocity  field  on  the  assumed  rigid-body  model 
must  be  taken  into  account.  At  this  point,  various  approximations  are  made 
to  simplify  the  modeling  problem.  One  such  approximation  is  given  in  [30] 
where  the  space  distribution  of  the  velocity  field  is  lumped  at  the  various 
stations  on  the  body.  The  resulting  equivalent  velocities  are  used  for  com¬ 
puting  aerodynamic  forces  and  moments.  In  [30],  the  penetration  effects  of 
the  wind  gusts  are  also  considered.  The  treatment  with  such  depth  produces 
a  relatively  complex  model  for  the  air-mass  velocities;  it  is  recommended 
for  cases  where  the  body  velocities  with  respect  to  ground  are  relatively  low 
and  the  bx  dy-bending  modes  are  significant. 


In  ADAPS,  a  simpler  model  canbe  used  because  of  relatively  high  aircraft 
and  weapon  velocities  and  rigid-body  model.  It  is  assume  a  that  the  velocity 
field  is  linearly  distributed  about  the  eg  of  a  rigid  body.  The  overall  influence 
of  .he  motion  of  air  mass  is  accounted  for  by  Va  and  its  space  gradient  matrix 
Bva 

d  evaluated  at  the  eg  of  a  rigid  body.  It  can  be  shown  [16,  31]  that  the 

effects  of  space  gradients  can  be  conveniently  taken  into  account  by  use  of 
the  equivalent  angular  velocity  vector  defined  by 


UJ  =  UU  +  U) 

a  w 


(4.  35) 


where  w  is  the  synthetic  angular  velocity  vector  corresponding  to  the  gradient 
effects  of  the  moving  air  mass. 

In  the  following,  a  simplified  model  for  the  equivalent  linear  and  angular 
velocities  is  presented. 


Modeling  for  the  Equivalent  Linear  and  Angular  Velocities 


First  the  modeling  for  the  mean  wind  is  presented.  Then  modeling  for  the 
turbulence  (i.  e. ,  gust)  is  treated. 


Mean  Wind  Model  --As  shown  in  Figure  16,  the  mean  wind  is  described  by  its 
magnitude  and  its  orientation  in  the  earth-fixed  axes  as  follows: 


(4.  36) 


Figure  16.  Description  of  the  Mean  Wind 
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where 


Then  along  the  body  axes  the  mean  wind  has  the  following  components: 

w  =  E  wg  (4.40) 

where  E  is  the  transformation  matrix  from  earth  to  body  axes  as  defined  by 
equation  (3. 14)  in  Section  III. 

Tlie  magnitude  V  of  the  mean  wind  is  assumed  to  be  altitude- dependent 
according  to  the  following  simplified  functional  relationship 


(4.  41) 
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where 


h  =  the  altitude  of  interest 

h  =  a  reference  altitude  at  which  the  mean  wind  speed  is  known 

o 

V  =  the  mean  wind  speed  at  reference  altitude  h 
o  o 

e  =  an  empirical  exponent  whiui  expresses  the  thermal  stability 
conditions  of  the  atmosphere  and  the  degree  of  roughness  of 
the  surface  beneath 

For  slightly  unstable  air,  typical  values  of  e  range  from  0. 12  for  smooth 
surfaces  (such  as  deserts)  to  0.  38  for  very  ro-gh  terrain.  In  ADAPS,  a 
constant  value  of  0. 25  is  used,  representing  average  conditions.  Equations 
(4.  39),  (4. 40)  and  (4. 41)  constitute  the  mean-wind  model  used  in  ADAPS. 

Gust  (Turbulence)  Model  --  The  gust  model  used  in  ADAPS  is  the  form 
attributed  to  Dryden  with  the  coefficients  specified  in  Ref.  31.  In  this 
form,  it  is  assumed  that  gust  velocities  are  locally  isotropic  (i.  e. ,  locally 
invariant  with  respect  to  position  and  orientation)  and  that  time  variations 
are  statistically  equivalent  to  distance  variations  in  traversing  the  gust  field. 

The  translational  gust  velocity  vector  is  defined  as 


(4.  42) 


The  power  spectral  densities  for  the  translational  gust  velocity  components 
are  given  by 
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where 


Q=  2IL 
a  X 

=  spatial  frequency,  (rad/ ft) 

(4.  44) 

X 

=  wavelength,  (ft) 

Li 

=  scales  (ft) 

0. 

1 

=  the  root  mean-square  gust  velocities  (ft/sec) 

i 

£ 

* 

It 

The  mean-square  gust  velocities  and  the  scales  are  related  to  each  other 
through  the  following  set  of  equations: 


o2  o2 

U  _  V 

L  “  L 

U  V 


w 


w 


(4.  45) 


The  quantities  appearing  above  have  the  following  altitude  dependence: 


100  <h<  1750  ft 


L, 

=  h 

w 

=  Ly  =  145. 0  h1^3 

Lu 

(4.  46) 

h  > 

1750  ft 

Lw 

=  L  =  L  =  1750  ft 
u  V 

=  5.  25  -  log- 


1.25 


w  —  ”°10  (10,000j 

For  0  <  h  <  100,  the  value  of  h  =  100  is  used  in  the  above  equations. 
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The  gradient  effects  of  gust  velocities  are  considered  as  explained 
earlier  by  defining  angular  gust  velocities  as  follows: 

/  3w  \ 

/n  V  / - g  \ 


- g 

sy 


(4.  47) 


The  spectra  for  the  angular  gust  velocity  vector  are  given  by 

r  .i/sT 


0  <0> 
Pg 


l)  <°-8>hnr 

Lw/  Ci+jfn)2] 


§  (fi)  =  0  (Cl) 

g  Mg 


- . ;  -  0W  (n) 

[i+fo2]  Wg 


(4.  48) 
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where  b  =  wing  span,  (ft) 

Random  velocities  with  above  spectra  are  obtained  j y  passing  a  gaussian 
random  "white"  noise  through  a  linear  system  with  a  proper  transfer  function 
c(s) 

^  G(s)  * 


0Q(w) 


It  is  known  that 


0O<«>  =  jG(ju)  J  0.(U»)  (4.49) 

where  the  bars  denote  the  magnitude  of  the  complex  variable. 

The  power  spectral  densities  given  above  are  ratios  of  polynomials  in 
where  w  is  the  temporal  frequency  given  by 


ui  =  V  Q  rad/ sec 
St 


(4.  50) 


can  be  spectrally  factored  out.  This  process  yields  the  proper  transfer 
functions  as  follows: 


‘k.-.  Ai.i  j  v  n.  l. 
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Figure  17  shows  the  structure  of  the  equivalent  linear  and  angular  gust 
velocity  generator  model.  Figure  18  shows  its  state  representation,  and 
Table  II  gives  the  filter  coefficients. 


The  pole  locations  of  roll,  pitch  and  yaw  gust  filters  for  aircraft  are 
inversely  proportional  to  wing  span,  b.  The  corresponding  term  for  a 
weapon  would  be  its  diameter.  When  weapon  diameter  is  used  in  place  of 
wing  span,  the  magnitudes  of  poles  of  roll  pitch,  and  yaw  filters  become 
excessively  large.  Numerical  integration  (i.  e. ,  non- real-time  simulation) 
of  these  extremely  fast  dynamics  requires  very  small  integration  step  size. 
For  small  weapons,  the  space  gradient  effects  of  wind  gust  (i.  e. ,  roll, 
pitch  and  yaw  filters)  are  small.  For  this  reason,  these  filters  are  omitted 
in  the  simulation  of  weapons  in  ADAPS. 

Steady-State  Output  Variances  of  the  Dryden  Gust  Filter  —  In  the  following, 
only  side  gust  filter  covariance  analysis  is  presented.  Others  follow  the 
same  pattern. 


! 

i 


i 


The  transfer  function  for  the  side  gust  is  given  as  [see  equation  (4.  51)] 


This  can  be  written  as: 


Gvg(s> 


s+b 


(s  +  a) 


(4.  51a) 


(4.  55) 


Equation  of  the  Gust  Generator 


Finally  using  (4.  57)  in  (4.  64)  results  in 


(4.  65) 


This  shows  that  in  order  to  obtain  variance  o  at  the  output,  the  input 
noise  r\  must  have  a  variance 


(4.  66) 


instead  of  unity,  or  the  gain  element  k  in  transfer  function  (4.  51a)  should 
be 


k  =  c 


yw 

vT  V 


(4.  67) 


with  a  unity  input  covariance. 

Subroutine  WINDK  computes  the  coefficients  of  the  gust  filter  having  the 
Dryden  spectrum  and  the  components  of  the  mean  wind  along  the  body  axes  as 
a  :  unction  of  altitude. 


SECTION  V 

DEVELOPMENT  OF  MEASUREMENT  SYSTEM  MODEL 


Sensors  on  board  an  aircraft  can  be  divided  into  two  groups:  instruments 
for  the  automatic  control  of  vehicle  motion  and  instruments  for  the  avionics 
(i.  e. ,  fire  control  navigation,  etc. ). 

In  general,  the  readings  (i.  e. ,  observations)  of  sensors  on  board  an 
aircraft  depend  upon  where  and  how  the  sensors  are  mounted  with  respect 
to  the  body  axes  of  the  aircraft  (geometry).  There  may  also  be  dynamics 
associated  with  sensors. 

To  provide  an  analysis  tool  by  which  the  effects  of  various  measurement 
points  and  methods  can  be  investigated,  both  aspects  are  considered  in  the 
development  of  a  measurement  model  in  ADAPS. 

In  the  following,  a  geometric,  model  for  the  overall  measurement  system, 
is  developed  first  (i.  e. ,  observation  equations).  Then  the  dynamics  of" 
sensors  are  treated.  Throughout  this  development  it  is  assumed  that  the 
body  on  which  instruments  are  mounted  is  rigid.  Aeroelastic  effects  for  the 
measurements  and  controls  are  beyond  the  scope  of  the  present  program. 


DEVELOPMENT  OF  A  GEOMETRIC  MODEL  FOR  MEASUREMENTS 

The  basic  vector  quantities  which  may  be  measured  are: 

•  Control  surface  deflection  vector: 

x6  =  col  <6a.  5g»  cr.  6sp;  (deg) 

•  Linear  acceleration  and  velocity  vectors: 

t  ,  o 

Xj  =  col  (u,  v,  w)  (ft/sec  ) 

Xj  =  col  (u,  v,  w)  (ft/sec) 

•  Angular  acceleration  and  velocity  vectors: 

x2  =  col  (p,  q,  r)  (rad/sec2) 
x2  =  col (p,  q,  r) (rad/sec) 

•  Angular  position  vector 

(i.  e. ,  body  attitude  with  respect  to  earth-fixed  axes) 

Xg  =  col  (0,  if  ,0) 
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•  Translational  position  vector 


x4  =  co1  (xe'ye,ze)  (ft) 

As  is  shown  in  the  subsequent  subsections  and  illustrated  in  Figure  20, 
instruments  sense,  in  general,  the  nonlinear  combinations  of  these 
quantities.  This  is  referred  to  as  the  geometry  of  the  measurements 
or  in  state  space  terminology,  the  observation  relations. 


MEASUREO 

VARIABLES 


Figure  20.  General  Structure  of  Perfect-Measurement  Model 


In  the  following,  first  the  geometry  of  the  measurements  for  the 
automatic  control  of  vehicles  motion  are  given.  Subsequently,  a  geometric 
model  for  air  data  measurements  is  presented  for  completeness.  Then 
a  simple  geometric  model  for  fire  control  measurements  (i.  e. ,  radar 
measurements)  is  developed. 


Geometry  of  Control  Measurements 

Attitude  Measurements  --  The  geometric  model  of  very  simple  forms  of  free 
gyros  (i.  e.,  vertical  and  directional  gyros)  are  given  in  [15 j.  In  many 
advanced  vehicles,  however,  more  complex  attitude  and  direction-sensing 
instrumentation  is  used.  Models  for  those  which  are  used  in  fire  control 
system  are  treated  later. 


js  ^ :v ti-r-  „ -  ^\<rr^Mk^!^'i*-r?v ^*^vr'*x*?*'#'.^w*?^fi 7fpv' 


^  7"  ^  7«f  'Vi'^11  •r^'1  VB'  v 
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The  relation  between  the  attitude  state  x„  and  the  sensed  (i.  e., 
measurable)  attitude  ygg  is  given  by  ° 

y3s  =  hx3(x3*  EtJ  EftX3  (5.1) 

where  is  a  fixed  transformation  matrix  which  takes  into  account  the 
nonaligmnents  effects  of  the  gyro  axes  and  is  normally  equal  to  an  identity 
matrix. 

Velocity  and  Acceleration  Measurements  -- 

Linear- Velocity  Measurements  --  The  geometric  model  of  linear- 
velocity  measurements  is  given  in  Figure  21.  In  this  model,  it  is  assumed^ 
that  the  origin,  0V,  of  the  instrument  axes  system  is  located  by  a  vector  Arv 
from  the  origin,  0,  of  the  body  axes.  It  is  further  assumed  that,  ijrv,  0V  and 
0V  are  fixed  Euler  angles  from  body  to  instrument  axes,  and  Ev  is  the 
corresponding  transformation  matrix  as  described  in  Section  III. 


Figure  21.  Axes  Systems  for  Linear  Velocity  and  Acceleration 
Measurements:  (a)  Body  Axes,  (b)  Instrument  Axis 


The  linear  (i.  e. ,  translational)  velocity  of  the  point  0V  in  body  axes  is 
given  [16]  by 

vSB  =v+53xArv  (5-2) 


where 

v  =  linear  velocity  of  eg  with  respect  to  earth  axes  in  body-axes 
system 

u5  =  angular  velocity  of  body  with  respect  to  earth  axes  in  body-  axes 
system 

^ry  =  position  vector  from  0  to  0y  in  body-axes  system 
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As  described  in  Section  III,  the  matrix  equivalent  of  (5.2)  is 

vSB  =  v  +  WArv  (5.3) 

where 


USB 

0  -r  q 

Ax 

V 

VSB  " 

VSB 

,  W  = 

r  0  -p 

,  Ar  = 

V 

Ayv 

-q  p  0 

_Azv_ 

Transforming  (5.3)  from  body  to  instrument  axes  system  yields  the 
measurable  velocity  components  along  the  instrument  axes 

vg  =  Ev  [v  +  W  Ary]  (5.  5) 


Equation  (5.5)  constitutes  the  geometric  model  of  the  linear  velocity  observa 
tions.  It  is  noted  that  the  measurable  velocity  vector  is  a  linear  combin¬ 
ation  of  the  linear  and  angular  velocities  in  body  axes;  Ey  and  Ar  are  the 
velocity  observation  parameters.  v 

In  state  vector  notation  (5.  5)  becomes: 


yls  =  \  (xl'  x2'  ArV  V 

where 


(5.6) 


hXl=EV[xl+LW(x2):l  *V} 

and  W(x2)  is  defined  by  (5.4). 


(5.7) 


Linear- Acceleration  Measurements  --  The  accelerometers  are 
assumed  to  be  located  and  oriented  in  much  the  same  way  as  the  velocity 
instruments.  In  the  following  the  geometric  relations  which  express  the 
accelerometer  observations  in  terms  of  the  vehicle  body  axes  acceleration" 
are  developed  similar  to  previous  section. 


Differentiating  (5.2)  in  a  rotating  frame  of  reference  (i.  e., 
yields  the  following  equation: 

vs  -  ^[v+uuxAra]+ujx  [v  +  3  x  AraJ 
Noting  that  ArQ  is  constant,  one  obtains 

cl 


v  =  -££+  JliLx  Ar_  +ujxv  +  ujx(5ux  Ar  ) 
s  ot  ot  a  a 


body  axes) 

(5.8) 

(5.9) 
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where 


ir  =  acceleration  of  point  0  where  accelerometer  is  located 
fa  expressed  in  body  axes' 

V  =  acceleration  of  mass  center,  0,  expressed  in  body  axes 
A  r  =  position  vector  from  0  t-  0_,  expressed  in  body  axes 

»  a 

The  matrix  equivalent  of  (5.  9)  is: 

v0  =  v  +  W  v  +  (W+W2)  Ar  (5.10) 

S  a 

Now  let  i)  ,  0  ,  and  0  be  the  fixed  Euler  angles  from  body  to  accelo- 
meter  axes,  Mid  Ea  be  the  corresponding  transformation  matrix,  then 
the  observable  acceleration  components  along  the  accelerometer  axes 
are  given  by 

a  =  E  [(a+W  Ar  )  +  (Wv+W2  Ar  )]  (5.11) 

o  a  a  a 

Equation  (5. 11)  constitutes  the  geometric  model  of  the  linear  accelerometer 
observations.  It  is  noted  that  the  observable  acceleration  vector,  a  , 
is  linearly  dependent  to  body  axes  accelerations,  but  also  contains  s 
nonlinear  combinations  of  the  velocities.  Ea  and  Ara  are  the  acceleration 
observation  parameters.  A 

In  state  vector  notation  (5. 11)  becomes 

yis  =  hx  x2*  xl'  x2'  Ara'  EaJ  12* 

where 

hx1  =  ^{*1  +  [W<i2)  ]  Ara  +  W(x2>xl  +  [W(*2>'2  ara}  (5. 13) 

Angular- Velocity  and  Acceleration  Measurements  --  Here  it  is 
tacitly  assumed  that  the  axes  of  the  instruments  which  measure  the 
angular  velocity  and  acceleration  with  respect  to  nonrotating  earth  are 
fixed  with  respect  to  the  body  axes.  It  should  be  noted  that,  in  radar-based 
measurements  which  are  treated  later,  the  radar  axes  system  on  whicn 
rate  sensors  are  mounted  moves  with  re~pect  to  the  body. 

Let  E  and  E^  be  fixed  transformation  matrices  from  body  to 
angular  rate  and  angular  accelerometer  axes  respectively.  Then  the 
observed  angular  rates  and  accelerations  are  given  by 

a.  =  E  w 

ti  U) 


Here  E^and  E^  are  the  observation  parameters.  In  state  notations  these 
equations  become 


where 


y2s  =  W 

H*  =  W 


(5.14) 


VCEJx2 

hj,  =  CE^  ]  lt2  (5.15) 


Air  Data  Measurements 

Three  measurements  are  made  on  the  air  data:  (1)  stu.lc  pressure  (2) 
total  pressure  and  (3)  total  temperature.  These  measure^.  *  '  are  used  by 

the  air  data  computer  to  obtain,  among  other  things,  the  aV  :  de,  altitude 
-ate,  mach  number  and  airspeed. 

In  the  following,  the  geometry  of  the  air  data  measurements  are 
developed  via  the  physics  of  the  variables  which  are  observed. 

Static  Pressure  Measurement  —  The  observed  static  pressure  is  a  non¬ 
linear  function  of  the  altitude.  The  static  pressure-altitude  relation  is  de¬ 
rived  from  the  barometric  equation  which  may  be  expressed  in  the  following 
form  [38]* 


g  W 

d  loEe  Ps  =  "  RT  ^ 
and 

T  =  T  +  t 
o 

where 

h  =  altitude 

g  =  acceleration  of  gravity 

W  =  molecular  weight  of  air 
m 

R  =  gas  constant 

T  =  absolute  temperature  (Kelvin) 

t  =  temperature  (centigrade) 

T  =  constant 

o 

Approximately,  one  can  write  from  (5.  16) 


(5.  16) 
(5. 16a) 
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or  in  state  space  notation 


Ps  =  Pb  e 


kl  ' 

T~ C  X4 


<5. 18) 


where 

c  =  col  (0  0  1) 

and  pb<  kj  are  the  observation  parameters. 

Total  Pressure  Measurement  -  -  The  observed  total  pressure  is  connected 
to  the  static  pressure  and  the  speed  of  air  with  respect  to  body,  through  the 
following  relationship: 

r 

PT  =  Ps  1  +IF  T  (5.19) 

a 


(5. 19) 


where  y  =  1.  4 


and  V  =  airspeed  with  respect  to  body 
a  =  speed  of  sound 

At  this  point,  it  is  convenient  to  introduce  the  Mach  number  parameter 
defined  by 


(5.20) 


Temperature  observations  are  functions  of  this  parameter. 

The  speed  of  sound  is  related  to  the  temperature  as  follows: 


a=V\Vv* 


(5.21) 


r  <a0 

where  k  *  y~ ^ 
o  p  f 
ro  o 


(5.22) 


Nov,'  combining  (5.21)  with  5.19)  and  noting  that 


Va  =  V  -  W 

o 


(5.23) 
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One  can  write  in  terms  of  states,  the  following  observation  equation: 


PT  =  Ps 


1  + 


ktfrj^ia*'  <xla> 


(3.5) 


(5.24) 


where 


x.  =  (x  -  x  ) 
la  1  iw 

k.  =  k  T 
1  o  o 


(5.25) 


Total  Temperature  Measurement  --  The  total  temperature  is  related  to 
the  static  temperature  t,  by  £38] 

tg  =  (1  -  .2M2)  (1  -  .004M2)t  (5 

and  is  observed  by  a  resistive  sensor  obeying  the  Callendar-Van  Dusen 
equation 


RT  =  g(ts)  (5.27) 

where  the  function  g  is  a  second-degree  polynomial  in  t  .  Thus  in  (5.24)  the 
parameter  t  (i.  e.,  temperature)  is  indirectly  observed  find  this  observation 
is  described  by  (5.26)  and  (5.27). 

In  summary,  the  air  data  observation  vector  is  connected  to  the  states 
as 


pss  =  hps(x4*  pb'  0 
pTs  =hpT(V  xlw-PSS'  0 

**Ts  hTT(xl’  x  lw*  ^ 


(5.26) 
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where 


/' 


ps 


p.  e(T  +t)  C*  X4 

D  O 


(h  => 
ps 

<Vxiw),(VxuJ 

h  rr,  = 

pT 

‘  *  *0<Vl> 

II 

JC, 

g<U 

-,3.5 


(5.29) 


Fire  Control  Measurements 

Description  of  Fire  Control  Measurement  Model  --  Fire  ccntrol  measure- 
ments  involve  body- mounted  sensors  also.  ?rom  the  readings  of  these  sen¬ 
sors,  the  target's  relative  position  and  velocity  with  respect  to  aircraft 
are  derived  [39], 

The  radar  measurement  model  considered  here,  is  illustrated  in 
Figure  22.  In  this  model  the  relative  position  of  target  is  defined  by  the 
vector 


RP 


R 


R 


(5.30) 


whore 

R 


R 


magnitude  of  the  position  vector  of  target  relative  to  point  Or  on 
aircraft 


i)ir  =  azimuth  angle  of  position  vector  with  respect  to  body  axes 
0R  =  elevation  angle  of  position  vector  with  respect  to  body  axes 

It  is  assumed  that  the  relative  position  vector  yj{p>as  defined  above  is  observed 

using  a  radar  device  which  is  located  by  a  vector  from  the  origin,  O,  of 
the  body  axes.  The  rad^r  axes  are  ass  uned  to  be  oriented  by  synchros  so 
that  the  antenna  (i.c.  ,  vector)  always  points  to  the  target.  It  is  further 
assumed  that  the  angular  velocities  qp  and  rR  of  radar  axes  with  respect  to 
earth  are  observed  by  antenna- mountea  rate  gyros,  and  the  rate  of  change  of 
Rpj  is  observed  by  doppler  shift.  This  describes  the  fire  control  measure¬ 
ment  model  used  in  ADAPS.  In  the  following,  the  geometry  of  the  radar- 
based  measurements  are  developed  parallel  to  the  previous  subsections. 
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Since 


drR  _  - 

“3T  "  “vr 


(5.37) 


where  ^  is  the  linear  velocity  of  Op  with  respect  to  earth-fixed  frame 
origin  0T.  One  obtains  from  (5.36)1ind  (5.37): 


-*  6r„  — •  -* 

*VR  =  -R  +  <*R  x  r 

fit 


In  matrix  notation 


fir 


where 


-vR  =  6tS+WR'R 


(5.38) 


(5.39) 


[R 

rR=|° 


fir 

It 


R 


|R 

=  I  0 


Using  (5.  40)  in  (5.  39)  yields 
•  •  “* 

-R 
-R  r 

R, 


w  = 
R 


0  *  ~rR  qR 


R. 

l  qR ; 


0 

PR 


R 


(5.40) 


(5.41) 


It  is  noted  that  R,  R,  and  angular  rates  qR  and  r  of  radar  axes  with 
respect  to  earth  axes  are  the  observed  quantities/ 

Now,  the  linear  velocity  of  Or  can  be  expressed  in  terms  of  the  velocity  of 
eg  in  body  coordinates:  K 


VRB  =  v+U)  X  ArR  (5.42) 

Let  Er  =  Er(0„,  ♦  _)  be  the  transformation  matrix  from  body  axes  to  radar 
axes.'lE  hroOT&in&d  from  equation  (3. 14)  of  Section  III,  by  letting  0  =  6P, 
i|r  =  fR  arflJ  <t>  =  o).  K 

Then  in  matrix  notation  the  linear  velocity  of  Or  in  radar  axes  becomes 

vR  ~  £ER<eR*  Yr)  3  Cv  "■  ArR]  (5.43) 

where  vR  is  given  by  (5.  41). 
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Defining  the  velocity  o! 

Fr 

A 


bservations  by 


yRV 


R 


(5.44) 


one  can  write  us-ng  (5.41) 
yRV  ■  D(R)  vr 


where 


D(R)  = 


-10  0 
1 
R 

0  0 


0 


(5.45) 


(5.46) 


Finally  substituting  (5.43)  into  (5.45)  yields  the  observation  equations  in 
the  form 


where 


yRV  hRV  (xl’  X2'  yRP'  ArR) 

hRV  =  D(R)  [Er(0r,  xl/n)]  [x1  +  W(x2)  ArR  1 


(5.47) 

(5.48) 


This  finishes  the  treatment  of  the  observation  equations  of  the  measure¬ 
ment  system  considered  in  ADAPS. 

In  the  subsection  that  follows,  the  development  of  a  model  for  the 
dynamics  of  sensors  which  read  the  above  observations  is  briefly  presented. 


DEVELOPMENT  OF  A  DYNAMICAL  MODEL  FOR  SENSORS 

Almost  invariably,  the  dynamics  associated  with  each  measured -scalar 
signal  is  of  second  order.  Therefore,  the  overall  system  order  increases 
very  rapidly  when  the  number  of  measured  signals  increases.  To  overcome 
this  difficulty,  the  sensor  dynamics  with  poles  lying  outside  of  the  signifi¬ 
cance  circle  of  radius  R  on  the  complex  plane  as  shown  in  Figure  23  are 
ignored  in  the  dynamicar representation  of  the  overall  system.  However, 
their  positions  are  checked  after  the  optimal  gain  loop  is  closed,  because  of 
the  sensitivity  of  high-frequency  open-loop  poles  to  feedback. 

•  ttl 

Figure  24  is  the  block  diagram  of  the  i  sensor  dynamics. 
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Figure  23.  '  Dynamical  Significance  Radius 


yio 


In  this  diagram,  y.  and  y.  correspond  to  the  i*  scalar  observable  and 
mea^red  signals  i^spectiWily;  x^  and  x2.  are  the  state  components  of 
the  i  sensor. 

The  dynamics  of  each  sensor  are  identified  by  five  coefficients,  a.., 
a„.,  b..,  b„.,  and  d.,  where  i  is  the  sensor  index.  These  coefficients  11 
cotresfyo  na  *to  the  oiltput-frobenius  implementation  of  the  sensor  transfer 
function.  For  those  cases  in  which  the  transfer  coefficients  (i.  e. ,  di)  are 
zero,  this  representation  provides  output  as  the  first  component  of  sensor 
state.  This  finishes  the  development  of  a  dynamical  model  for  sensors. 

Figure  25  shows  in  detail  the  kinematics  and  dynamics  of  measurements. 
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SECTION  VI 


LINEARIZATION  OF  THE 
WEAPON  DELIVERY  PROCESS 


The  weapon  delivery  process  is  a  nonlinear  stochastic  phenomenon.  To 
analyze  the  nonlinear  equations,  practically,  they  are  linearized  about  the 
nominal  path.  The  nominal  path  considered  here  corresponds  to  a  dive- toss 
maneuver  and  consists  of  essentially  three  phases:  (a)  dive,  (b )  pull-up,  and 
(c)  free-fall,  as  shown  in  Figure  26. 

•V 

The  development  of  a  model  for  the  release  transient  phase  is  outside  the 
scope  of  this  work.  This  part  of  the  nominal  trajectory  is  taken  into  account 
in  a  simplified  manner,  by  introducing  an  independent,  additive,  stochastic 
error  on  the  initial  condition  for  the  free-fall  trajectory. 


Figure  26,  Nominal  Trajectory  for  Dive-Toss  Maneuver: 

(a)  Dive  Phase,  (b)  Pull-up  Phase,  (c)  Free- 
Fall  Phase,  (d)  Release  Transient  Phase 
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LINEARIZATION  PROCESS 


It  is  assumed  that  the  dive  and  pull-up  maneuvers  are  carried  out  with  a 
fixed  thrust  level,  and  by  controlling  the  elevator  deflection.  The  missing 
states  and  parameters  along  the  paths  (a)  and  (b)  are  obtained  either  by 
solving  a  set  of  trim  equations  developed  for  the  case  at  hand  or  by  a  soft 
flight-path  controller.  The  trajectory  (c)  is  developed  by  integrating  the 
six  degree -of-freedom  free-fall  equations  of  an  iron  bomb  ' 

The  linearization  process  along  the  paths  (a)  and  (b)  yields  the  state 
transition  and  the  input  matrix  pair  (F,  G)  of  the  aircraft.  Which  are  used  to 
design  an  optimal  weapon  delivery  controller.  The  linearization  along  the 
path  (c)  yields  the  sensitivity  matrices  which  are  used  in  translating  the 
dispersion  errors  at  bomb  release  to  the  Impact  of  the  bomb  on  the  target 
as  well  as  errors  which  occur  during  the  free-fall. 

The  numerical  development  of  the  ideal  nominal  trajectory  is  illustrated 
in  Figure  27  for  the  algebraic  trimming.  The  process  starts  at  t  »  to  with 
altitude  h0  and  range  x0.  First  a  trimming  is  performed  as  outlined  in 
Section  VII;  then  a  linearization  is  performed.  This  sequence  continues 
with  every  AT  seconds  along  the  dive  path  until  the  pull-up  altitude,  hpu,  is 
reached.  Similar  steps  are  carried  out  along  the  pull-up  trajectory,  until 
the  release  time,  tr,  is  reached.  Subsequently,  integration  of  the  six 
degree-of-freedom  equations  of  motion  of  the  bomb  is  carried  out  N  steps 
with  At  ■  AT/N  time  interval.  Then  a  linearization  is  performed  at  t  * 
tr  +  AT.  The  process  Is  continued  until  the  impact  plane  is  reached. 

The  equations  which  describe  the  general  motion  of  a  rigid  body  are 
given  in  Section  III.  Many  problems  of  rigid-body  motion  involve  only  small 
disturbances  (i.  e. ,  perturbations)  from  steady  or  quasi-steady  flight  condi¬ 
tions.  In  the  following,  the  assumption  of  small  disturbances  from  reference 
flight  conditions  is  used  to  reduce  the  equations  from  nonlinear  to  linear  form. 


DEVELOPMENT  OF  THE  PERTURBATION  EQUATIONS 

In  the  state  vector  notation,  the  general  equations  of  motion  of  a  rigid 
body,  as  developed  in  Section  HI,  is  described  by  a  nonlinear  vector 
differential  equation  of  the  form 

x  =  f  Cx(t),  u]  (6. 1) 

where  f  is  real,  continuous  and  has  continuous  first  order  partial  deriva¬ 
tives  with  respect  to  xi,  i  =  1, 2, . . . ,  n;  and  uj,  l  =  1,  2, ... ,  r;  in  a  region 
of  (  x,  u)  space  which  contains  the  solution  curve  (x(t),  u)  with  t0  *  t  <  tj. 
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Let  (x,  u)  and  (x,u)  be  two  neighboring  pair  of  solution  curves 
satisfying  16. 1).  Define 

?  =  x  -  x 

and  ~  _  (6. 2) 

T1  =  U  -  U 

Also  let  the  matrices  with  columns 

»  t  *■  I#  i  •  •  #  n, 

(x(t),u) 


be  denoted  by  F~  (  x(t),  u),  Fu  (  u(t),  u),  respectively.  Then  from  (6.  2)  it 
follows  that  L40J 


Yf  -  f  (  S+x(t),  h+u)  -  f  (  x(t),  u) 


(6.3) 


or 

df- =  Fx(  x,u)§  +  F^  (x,u)r)  +  o(|§j)  +  o(jn|)  (6.4) 

where  o(g)  is  a  vector  such  that 
lim  o (e)  _  n 

e-0  e  "  0  '6*5' 

When  the  last  two  terms  in  (6. 4)  are  omitted,  there  occurs  the  linear 
system 

=  Fx(  x(t),u)y  +  Fu  (  x(t),u)v  (6.6) 

with 

y(tQ)  =  x(tQ)  -  x(tQ)  (6.7) 

which _is  called  the  first  variation  of  (6. 1)  with  respect  to  the  solution 
(x(t),  u).  It  is  also  called  the  variational  equation  of  (6. 1).  The  first 
variation  determines  the  dependence  of  solutions  on  the  intial  conditions 
and  parameters.  It  also  determines  in  some  cases  the  nature  of  the 
stability  of  the  solutions  (x,  u)  of  (6. 1). 


vi*> i'tji r-*%*  *g 


In  engineering  practice,  the  procedure  described  above  is  called  the 
linearization  process,  and  the  variational  equation^  (6. 6),  is  called  the 
linearized  equation  of  motion.  Also  the  solution  (x,u)  is  referred  to  as  the 
nominal  solution  or  the  reference  trajectory.  It  follows  from  (6.  6)  and  the 
definitions  of  Fv_and  Fu  that,  in  order  to  carry  out  the  linearization: 

(a)  the  solution  (x,  u)  must  be  specified  on  an  interval  tQ  £  t  £  tf,  and 

(b)  the  first  partials  of  f(t,  x,  u)  with  respect  to  xi  and  uj  must  be  developed. 

In  the  following,  first  the  development  of  the  nominal  solution  (x,  u)  is 
given  then  the  development  of  first  partials  is  presented.  At  this  point  a 
few  words  on  the  notation  for  small  disturbances  'i.  e. ,  perturbations)  are 
in  order. 

Usually,  perturbations  of  velocity  and  orientation  variables  are  desig¬ 
nated  by  the  lower  case  symbols  for  these  quantities,  i.  e. , 


H 

fp\ 

(0 

/x  \ 
/  e  \ 

VJ 

»  1 

q  > 

6 

•  ye 

W'1 

ir  / 

l\|f  i 

ze 

Upper  case  symbols  are  used  with  a  subscript  zero  to  denote  the  reference 
values  of  these  variables.  Thus 


are  reference  values  for  linear  and  angular  velocity  components,  orienta¬ 
tion  angles,  and  positions.  Incremental  changes  in  aerodynamic  force  and 
moment  components  are  denoted  by  the  pertinent  symbol  with  a  prefix  A, 
e.g. ,  AX,  AZ,  AM,  etc. 


DEVELOPMENT  OF  THE  NOMINAL  SOLUTION 

The  way  the  nominal  solution  (x,  u)  is  developed  depends  upon  whether 
or  not  the  parameter  vector  u  is  controlled  or  uncontrolled.  The  term 
M  controlled"  here_implies  the  description  of  u  over  a  time  interval  sucl 
that  the  solution  (x,  u)  behaves  as  specified.  Uncontrolled  u  on  the  other 
hand  implies  generally,  disturbance  parameters  effecting  the  evolution  of 
the  solution. 

The  nominal  solutions  (i.  e. ,  reference  trajectories)  are  developed  by 
means  of  a  set  of  specified  reference-flight  conditions.  Reference -flight 
conditions  are  divided  into  two  groups:  free  reference-f Light  conditions 
and  controlled  reference -flight  conditions. 
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Free  re  fere  nee- flight  conditions  constitute  the  specification  of  an 
initial  value  of  the  state,  x(t0)  a  Xq  and  the  specification  of  the  free  parameter 
u  over  an  interval  tQ  ^  t  ^tf.  The  nominal  solution  in  this  case  is  developed 
by  integrating  (6. 1)  over  the  specified  time  interval, 

A  nominal  motion  of  an  iron  bomb,  released  from  an  aircraft  falling 
freely  under  the  influence  of  the  gravity  and  winds  is  an  example  of  this  case. 

The  set  of  controlled  reference -flight  conditions  consists  of  steady 
flight  conditions,  quasisteady  flight  conditions,  straight  flight  conditions,  and 
symmetric  flight  conditions  115]: 


Steady  Flight  Conditions  imply  a  motion  with  zero  linear  and 
angular  accelerations.  That  is 


dt 


U\ 


,w/ 


=  0,  and 


dt 


(6.8) 


Steady  sideslip,  level  turns,  and  helical  turns  are  examples  to  this 
kind  of  motion. 


•  Quasi  steady  Flight  Conditions  imply  a  motion  with  some  nonzero 
linear  and/or  angular  accelerations.  A  steady  pitching  motion 
(i.  e. ,  steep  dive)  which  is  described  by  the  quasisteady  flight 
conditions  ' 


is  an  example  to  this  case. 

•  Straight  Flight  Conditions  imply  a  motion  with  zero  angular  velocity 
components?  That  is 


r\ 

q  ■  0  (6.  10) 

\  r  / 

Steady  sideslips  and  dives  or  climbs  without  longitudinal  acceleration 
are  examples  of  this  kind  of  motion. 
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Symmetric  Flight  Conditions  imply  a  motion  in  which  the  body  plane 
of  symmetry  (i.  e. ,  xz  plane),  remains  fixed  in  space  throughout  the 
flight.  That  is  the  asymmetric  variables  are  all  zero: 


Wings  level  dives,  climbs,  and  pull-ups  with  no  sideslips  are 
examples  of  this  kind  of  motion. 

The  set  of  mathematical  consequences  associated  with  all  specified 
reference  flight  conditions  is  used  to  construct  the  nominal  trajectory. 

At  this  point  it  should  be  noted  that  some  of  the  reference  flight  conditions 
are  on  the  derivatives  of  the  state  variables  rather  than  on  the  state 
variables  themselves. 

The  process  o f  finding  the  values  of  involved  state  variables  and  inputs 
so  that  the  conditi  >ns  as  specified  by  the  equations  given  above  are  satisfied 
is  called  "trimming".  (The  process  of  trimming  for  a  particular  set  of  flight 
conditions,  i.e. ,  dive  and  pull-up  maneuver,  is  described  in  Section  VII, ) 


DEVELOPMENT  OF  THE  FIRST  PARTIALS 

Consider  a  general  rigid-body  motion  characterized  by 
x  =  f(x,  u,  w)  (6. 12) 

where 

x  =  state  vector  of  the  motion 
u  =  control  vector  of  the  motion 
w  =  disturbance  inputs  of  the  motion 

Equation  (6. 12)  can  be  decomposed  into  the  following  form  by  using  the 


Equations  (3.62),  (3.63),  (3.24)  and  (3.72)  of  Section  III: 

ij  =  -W(x2)x1+E(x3)ge+ i  fa<*1.*2.Xj.8,'*r>+£BTyT  (6.13) 

Xj  ■  -J‘1W(x2)J  Xj+ j"1  ma(x1,x2.x3,6.w)+ J’1  BTyT  (6.14) 

*3  =  C(x3)x2  (8.15) 

i4  =  E'(x3)Xj  (6.16) 
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where  the  subvectors  are  defined  as  follows: 


linear  velocity  vector 


angular  velocity  vector 


angular  position' vector 


w  = 


translational  position  vector  (flight  path  state) 


control  surface  deflection  vector 


effective  thrust  input 


wind  velocity  vector 


and  finally  fa  and  ma  are  the  aerodynamic  force  and  moment  vectors, 
respectively,  expressed  in  the  body  coordinates.  These  two  vector  functions, 
in  general,  do  not  have  analytic  form.  Their  values  as  functions  of  their 
arguments  are  supplied  in  the  form  of  tables.  In  this  work,  the  dependencies 
of  fa  and  ma  to  the  derivative  of  the  state  vectors  as  defined  above  are  ignored. 

If  a  quaternion  is  used  to  describe  the  angular  position  coordinates,  then 
equations  (6. 13)  through  (6. 16)  basically  remain  the  same  except  the  angular 
position  vector  becomes 


r*i 

x„ 


V 


3. 

V 


and  the  differential  equation  of  angular  position  can  be  written  as  [see 
Section  III,  Equations  (3,54)  and  (3.68)J 


x3  =  G(x2)  x3  =  G(x3)x2 


(6. 15a) 


where  G(xo)and  G(x3)  are  linear  functions  of  x3  and  X3  respectively.  The  state 
diagram  of  the  nonlinear  equations  of  motion  is  illustrated  in  Figure  28. 

Equations  (6. 13)  and  (6. 14)  describe  the  evolution  of  the  linear  and 
angular  velocity  vectors,  respectively.  Equations  (6. 15)  and  (6. 16)  describe 
the  evolution  of  the  angular  and  translational  position  vectors,  respectively. 
Two  approaches  are  available  for  the  development  of  partialsof  the  right-hand 
sides  of  (6. 13)  through  (6. 16):  (a)  mixed  (i.  e. ,  analytical  and  numerical) 
partial  differentiation  approach  or  (b)  pure  numerical  partial  differentiation 
approach.  Depending  upon  what  approach  is  used  for  obtaining  the  partials, 
linearization  will  be  referred  to  as  (a)  mixed  linearization,  and  (b)  pure 
numerical  linearization.  Both  approaches  are  presented  in  the  following, 
but  only  the  pure  numerical  linearization  approach  is  utilized  for  the  lin¬ 
earization  in  ADAPS. 

Mixed  Linearization  Approach 

It  is  intuitively  obvious  that  mixed  partial  differentiation  approach  for  the 
linearization  is  more  accurate  than  that  of  pure  numerical  partial 
differentiation  approach.  However,  as  will  be  seen  in  the  following  develop¬ 
ment,  it  requires  more  programming  effort  and  computer  memory  for  its 
implementation.  Observation  of  the  right-hand  sides  of  equations  (6. 13) 
through  (6. 16)  reveals  that  the  majority  of  the  analytical  terms  are  in  the 
form  of 


f(x.,  Xj)  =  F(x.)  Xj 


(6. 17) 


In  the  following  the  incremental  change  in  f(xi,  xj)  is  developed  for  various 
values  of  i  and  j: 

•  Case  1;  bee  Equation  (6. 13)]  Clearly,  an  incremental  change  in  f 
about  the  nominal  states  (x.,  x_)  in  terms  of  perturbations  fix.  and 
5x„  is  given  by  1 

v 

Sfix^Xg)  =  6  IW(x2)x13  =[W(x2)]|6x1  +  [SW(x2)>:i|  (6.18) 

A  A 


96 


On  the  other  hand 


*r  \xt  xi  aw .  ^  aw 

6L  Wx2)J  =  |fip  +  ^“|  5q  +  "§?  |  6r 


(6. 19) 


So  that 


aflXj,  x2)  =  gtwtxgjxj]  =  [W(x2)]  o  6  Xj  + 


(6.  20) 


Since  W (xg)  is  a  linear  function  of  xo,  indicated  partials  are  cc 
matrices.  When  the  components  of W(x2>  are  used  in  (6. 20)  in 
accordance  with  equation  (3.21)  of  Section  III,  one  obtains 


constant 


6[W(x2)x11  =  W(Xg)  5Xj  -  W(x2)  5x2 


(6.21) 


where 


0  -w  v 


W(xx)  =  w  0  -u  ,  W(x2)  =  r  0  -P  (6.21a) 

-v  u  0  -q  P  0 


0  -r  q 


The  same  result  can  also  be  obtained  from  the  vector  notation 
representation  of  f  by  noting  that 

?  =  wxv  =  -  vxw  (6. 22) 

so  that 

gf  =  w  |  x  gv  -  v  |  x  gw  (6.  23) 

o  o 

Equation  (6.  21)  is  matrix  representation  of  equation  (6.  23). 

Case  2: 

f  =  [j"V(x2)j]x2  (6.24) 


Following  the  similar  steps 


6f(x2,  x2)  =  6[J_1W(x2)J  x2J=  [j_1W(x2)jJo  6x 
+  (j 

•  Case  3: 


,X2)  * 


f(x2,  Xg)  =  Cg(x3)Jx2 
so  that 

6f(x2,  x3)  =  [G{x3)J  !  6x2  +  [6g(x3^Jx2| 


(6.25) 


(6.  26) 


(6.27) 


o 


But 


6(G(x3))  =  69  +  60  +  W  61)1 


(6.  28) 


So  that 


6f(x2,  Xg)  =  [G(Xg)-^  6x2  +  (^-x2|  ^  x2  ly^Xg  6xg  (6.29) 


•  Case  4: 


(6.  30) 


((Xj.Xg)  *  [E^Xgjlfj 

so  that 

6f(xr  x3)  =  E'(x3)]6  Xj  +  Jlf*"  Xj  xx  |  ■ xx  J  j>x3  (6.  31) 

In  summary,  the  linearization  of  the  analytic  terms  in  (6. 13)  through  (6. 16) 
yields  the  following  set  of  equations: 

6  Xj  =-[W(x2)]o6xj  -  (tfxjlfxjl^x,,  jlfs.lf  gel|f  Be)6xS 

,  ,  (6.32) 

+  SBT5yT+S  Afa<xrx2’x3-5’w) 


E  ~  77^7^77^77.  *J 


*• 

Si2  =  -  J''w(x2)J+  L'1  —■  Jx2|j 


"f^l^f  *,)  «*2 

'o 


- 1  A  -  1 

+  J  BT5yT  +  J  Ama(x1Jx2,x3,5,w) 


(6.  33) 


6xg  =  Cg(x3)Jo6x2  +  air  x2  I  90  x2  lyf  x2  j  &x3 

o 

&*4  =  Ce  (XgjJ^Xj  +  |-g-g*  x!  I  30"  x!  I  a^T  xi|o  5X 


(6.  34) 


(6.  35) 


where  the  matrices  W,  E,  and  G  are  defined  by  equations  (3.21),  (3.14)  and 
(3.  24)  of  Section  III,  respectively.  The  set  of  equations  (6.  32)  through  (6.  35) 
can  be  written  as 


mBT| 


[mAfa(xl’x2’x3’6’  w) 


-1 A  -1 

6  x2  J  Bj,  AyT-+  J  Am&(Xj,  x2,  Xg,  6,  w  ) 
+ 


(6.  36) 


This  finishes  the  linearization  of  the  analytical  terms.  What  remains  to  be 
done  to  finish  the  problem  is  to  compute  A  f  and  Am  appearing  in  (6.  36).  Since 
these  vector  functions  have  tabular  representations,  one  must  resort  to  a 
numerical  procedure  to  compute  them  in  terms  of  the  increments  in  their 
arguments.  At  this  point  there  exist  two  methods  of  approach:  (a)  the 
equations  of  the  force  and  moment  vectors  presented  in  Section  IV  can  be 
utilized  to  further  the  analytical  differentiation  process/  this  procedure 
finally  yields  to  the  stability  derivatives  of  the  aerodynamic  force  and 
moment  system;  or  (b)  direct  numerical  partial  differentiation  of  f  and  m. 


This  brings  us  to  the  problem  of  numerical  partial  differentiation  of 
vectors  with  several  arguments.  Since  the  size  of  a  vector  and  its  arguments 
are  immaterial,  the  computational  procedure  presented  in  the  section  that 
follows  for  the  pure  numerical  partial  differentiation  approach  applies  here 
as  well. 


Pure  Numerical  Linearization  Approach 


As  can  be  seen  from  the  equations  (6,  32)  through  (6.  35),  the  mixed- 
linearization  approach  requires  more  effort  in  its  computer  implementation 


and  it  does  not  completely  eliminate  the  numerical  differentiation  process. 
For  this  reason,  sacrifice  is  made  in  the  accuracy  of  computations  to  grossly 
simplify  the  linearization  process.  In  the  following,  a  set  of  candidate 
numerical  algorithms  is  presented  and  the  one  which  is  implemented  in 
ADAPS  is  discussed. 

Numerical  Differentiation  Algorithms  --  Numerical  differentiation  algorithms 
can  be  developed  either  by  (a)  differentiating  the  Lagrange  interpolating 
polynomial  of  a  table  function  or  (b)  using  the  discrete  approximation  to  the 
derivative  operator. 

By  differentiating  three-point  Lagrangian  interpolation  formulas  and 
evaluating  the  results  at  tabular  points,  the  following  derivative  formula  is 
obtained  141 3: 


-~-f'"(S),  xo-h<S<xQ  +  h 


(6.  37) 


If  it  is  known  that  If1  M(x)|  <  M3  in  the  interval  (x0-h,  xo+h)  and  if  all 
given  data  were  exact,  the  maximum  possible  error  in  the  calculation  of 
f'(xo)  would  be 


Mgh2 

1^*3 1  max  6 


(6.  38) 


On  the  other  hand,  if  each  of  the  ordinates  involved  is  in  error  by  ±  e,  then 
the  magnitude  of  the  corresponding  error  in  the  calculation  of  FCxq)  could  be 
as  large  as 


(6.  39) 


whereas  a  reduction  of  the  truncation  error  E3  would  generally  require  a 
decrease  in  h,  a  small  value  of  h  would  lead  to  a  large  possible  round-off 
error  R3  and,  conversely,  a  reduction  in  |  R3  jmax  would  generally  correspond 
to  an  increase  in  |  E3|max. 

A  reasonable  procedure  consists  in  determining  the  interval  h  such  that 
the  predictable  upper  bounds  on  the  two  errors  are  about  equal,  if  this  is 
feasible.  The  optimum  value  of  h  and  the  corresponding  maximum  total 
error  T*  are  then  found  to  be  [4l] 


h*  =  l,8e1/3M, 


T*  =  1. 1  s2'  3  Mg1/3 


(6.  40) 


By  using  a  discrete  approximation  to  the  derivative  operator,  the  derivative 
of  a  function  at  discrete  points  can  be  expressed  in  terms  of  the  differences  of 
the  function  at  tabular  points.  Let  s  be  a  differential  operator,  that  is 

sf(x)  =  f‘(x)  (6.41) 
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If  z  is  a  shift  operator  defined  by 


zf(x)  =  f(x+h)  (6. 42) 

Then  it  can  be  shown  by  using  Taylor  series  expansion  that 

z  =  ehs  (6.43) 

If  A  is  a  forward  difference  operator  defined  by  Af(x)  =  f(x+h)  -  f(x),  then  one 
can  write  z  =  ( 1+A )  and 

s  .  !^LJ.logn±A).|(A.i  42+>>>)  (6-44) 

In  terms  of  backward  differences,  z  =  (l-7)~* 

s  =  '^(l-v)  =  h<7+  \  v2+-**)  (6.45) 

If  the  expansion  is  truncated  after  the  first  term  one  obtains 

^«-J-(A)  (6.46) 

so  that 


Derivatives  using  more  data  points  can  be  obtained  from  (6.  44)  or  (6.  45). 

At  this  point  a  remark  is  made  about  algorithm  selection.  In  aircraft  or 
weapon  linearization  problems,  the  listed  increments  in  the  function  arguments 
are  much  larger  than  that  of  increments  that  can  occur  about  a  nominal 
argument  set  due  to  small  perturbations.  Consider  a  table  function 

f  =  f(li) 


where 

h(x,  u) 

|i  “  M(x,  u) 

_a\x,  u)_ 

is  the  data  argument  set  of  the  table  function.  Let 


b\x  = 


Ah 

AM 

Aa 
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be  the  set  of  increments  on  which  data  is  given.  Then  there  is  a  perturbation 
(§  x,  6u)  such  that 


|8U|  <  |4M|*|C1||6*|+  |C2||6u| 


(6.  48) 


This  fact  is  illustrated  in  Figure  29. 

h“ 


On  the  other  hand,  for  every 
represented  by 


M 

a 


in  the  i,  j,  k  data  cube,  £(h,  M,  a)  is 


f(h,M,a) 


l  +  <t1h  +  40M  + 

Ola  <5 


(6.  49) 


+  Jt12 h  M  +  ^23M  a  +  i31'a  h 


4>»  «V<  <%/ 


f  *12  3h  lVI  a 


where  i*s  are  functions  of  data  point  indices  i,  j,  k  only,  and  h,  M,  a  are 
linear  functions  of  (h,  M,  a).  This  representation  follows  from  the  multi¬ 
dimensional  linear  interpolation  algorithm  used  in  ADAPS  for  generating  a 
value  for  a  table  function  inside  the  interval  of  tabular  points.  It  follows 
from  (6. 49)  that  the  partials 

jaf  8f  9f| 

(TE'  m'  aa) 

o 


depend  only  upon  (a)  the  data  cube  in  which  the  nominal  parameter  vector  is 
located,  and  (b)  the  value  of  the  nominal  parameter  vector.  It  does  not 
depend  upon  the  size  of  the  perturbations  as  long  as  the  perturbation  cube  is 


Figure  29.  "Data  Increment"  and  "Parameter 
Perturbation"  Cubes 


inside  the  data  increment  cube.  Therefore,  we  may  conclude  that  a  two- 
point  linear  interpolation  algorithm  is  used  for  representing  a  tabular  function 
inside  the  data  interval,  and  when  small  perturbations  are  assumed  then,  the 
use  of  elaborate  differentiation  algorithms  is  not  justified  for  the  table  functions 
and  partial  differentiation  algorithms  given  by 

(6.  50) 


(6.  51) 


(6.  52) 

are  all  equivalent  and  moreover  independent  of  the  size  of  5«  provided  that 
1 6a |  <  Aa. 

For  pure  numerical  linearization  process,  however,  the  use  of  the 
central  differences  given  by  equation  (6.  52)  provide  better  accuracy  in  the 
linearization  of  analytic  functions  appearing  in  the  equations  of  motion.  For 
this  reason  (6.  52)  is  utilized  with  a  set  of  fixed  perturbation  increments,  in 
ADAPS. 

The  analysis  presented  up  to  this  point  is  implemented  as  the  subroutine 
LINK.  It  is  used  to  linearize  numerically,  the  nonlinear  state  equations  of 
motion. 


TRANSFORMATION  OF  THE  PERTURBATION  STATES 

A  few  words  are  in  order  here  abouc  the  state  component  assignment  to 
various  physical  quantities  in  the  dynamics.  The  standard  state  component 
assignment  for  the  nonlinear  equations  of  motion  has  been  defined  to  be 

x=  col(x1, 

=  col(u,  v,  w}p,  q,  r|  6,  0,  »  |  xg,  yg,  zg) 

The  resulting  state  perturbation  equation  of  motion  is  illustrated  in  Figure  30. 

When  dealing  with  the  linearized  equations,  the  state  components  may  be 
reordered  with  respect  to  longitudinal  and  lateral  dynamics  for  convenience. 
The  state  component  assignment  for  this  case  is  defined  as 

6x=  col(5x  ,  6h  ,  5u,  68,  6q,  6w  |  6  y,  6 Jr,  6r,  6v,  60,  6p)  (6.54) 

v  C 

The  resulting  state  perturbation  equation  of  motion  is  illustrated  in  Figure  31. 


X2*  X3^X4’  X5'  X6 1 X7'  V  X9*X10’  Xll*  X12* 
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The  gust  filter  states  (Section  IV)  are  also  grouped  with  respect  to 
longitudinal  and  lateral  ex  ,itations: 

w^cotfwj,  w?,  w3,  w5|w8,  w2,  v/4,  wg)  (6.57) 

The  resulting  state  equation  of  the  full-gust  filter  is  illustrated  in  Figure  32 
and  the  state  equatipn  of  the  reduced  gust  filter  (rolling  gust  terms  omitted) 
is  illustrated  in  Figure  33.  These  definitions,  produce  almost  upper-block 
triangular  transition  matrices.  They  are  obtained  by  permutational  similarity 
transformations.  This  operation  is  referred  to  as  "Shuffling”  of  the  linear 
data  and  is  carried  out  by  Subroutine  SHUFF. 

Besides  permutational  similarity  transformations,  as  explained  above, 
there  are  transformations  induced  by  the  various  selections  of  physical 
variables  as  state  components.  For  instance,  instead  of  selecting  v,  w) 
as  state  components,  one  may  choose  (V,  P,  a)  to  describe  the  velocity  vector. 


Consider  the  state  differential  equation  given  by  (6. 12).  Let  §  be  a  chosen 
state  vector  related  to  standard  state  vector  x  by  a  nonlinear  relation: 

?  =  g(x>  (6.58) 

Then 

?  =  Mx.u)  (6.59) 

where 

h(x’u)  =  r§x\  f(x’u)  (6.60) 

Let  us  assume  that  x,  2^  and  5  ere  evaluated  for  each  nominal  point  using 
(6. 12),  (6.  58)  and  (6.  59).  Then  the  partials 

(H!  ‘(Is)  '(!fM!£)  ’  andllrl  can  be  computed  along  the 
o  o  o  o  o 

nominal  trajectory  during  the  linearization  process  with  the  standard  state, 
as  described  previously. 


(6.61) 

(6.62) 
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Figure  32.  Full -Gust  Filter  State  Equations  with  Channeled 
State  Component  Assignment 


SECTION  VII 


DEVELOPMENT  OF  THE  NOMINAL  STATE 
AND  PARAMETER  (TRIMMING) 


The  way  the  nominal  state  and  parameters  of  a  rigid  body  in  flight  are 
developed  depends  on  whether  the  body  is  controlled  or  uncontrolled.  When 
a  prescribed  nominal  trajectory  to  be  generated  by  a  controlled  body  (i.  e. , 
aircraft,  guided  missile)  is  under  consideration,  one  can  find  it  either:  (a) 
by  simulation  with  an  autopilot  for  given  flight  conditions  and  reading  out  the 
nominal  values  of  state  and  control  parameters  during  the  flight;  or  (b)  by 
assuming  that  a  body  is  forced  to  fly  in  accordance  with  given. flight  condi¬ 
tions  and  computing  the  missing  nominal  states  and  parameters  which  produce 
that  flight.  The  first  technique  is  called  "trimming  with  an  autopilot"  and 
the  second  "algebraic  trimming. " 

When  a  body  is  uncontrolled  (i.  e. ,  iron  bomb,  bullet)  the  nominal 
trajectory  cannot  be  arbitrarily  specified.  It  is  obtained  by  integrating  the 
differential  equations  of  motion  starting  from  a  given  initial  condition. 

In  the  following  the  development  of  the  nominal  state  and  parameters 
along  a  prescribed  trajectory  are  treated  for  »  controlled  body  first  using 
the  algebraic  trimming  approach.  Then  trimming  with  an  autopilot  is  dis¬ 
cussed.  In  ADAPS  the  latter  approach  is  utilized. 


ALGEBRAIC  TRIMMING 
Consider  an  aircraft  represented  by 
x  =  f(x,  u) 
y  =  h(x) 

where 

x  =  state  vector 
u  =  control  vector 
y  *  response  vector 
The  response  rate  is  then  given  by 

y  =  H(x)  f(x,  u)  (7.3) 


(7.1) 

(7.2) 
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The  response  vector  may  consist  of  a  set  of  trajectory  variables  such  as 


+  V2  4-  W2 

\e-tan’1? 


(7.5) 


Most  often,  y  and  some  of  the  components  of  x,  and  x  are  specified  along 
the  nominal  trajectory. 


For  instance 


rvl  [a 

=  [rj  1° 


(7.  6) 


indicates  a  quasi- steady  motion  with  a  constant  acceleration  and  a  constant 
flight-path  angle  (steep-dive  bombing). 

The  trim  probl  m  is  to  find  the  missing  state  variables  and  controls 
such  that  the  error  defined  by 


e(t)  =  |H(x)  f(x,  u)  -  yd  | 


(7.7) 


is  small. 


Algebraic  Trimming  for  Dive  and  Pull-up  Traiecto 


The  two-phase  nominal  trajectory  considered  here  corresponds  to  a 
steady,  symmetric,  straight -path,  dive-bombing  maneuver  [7],  followed  by 
a  symmetric,  steady-pitch  (i.  e. ,  constant- g)  pull-up  maneuver. 

In  the  first  part  of  the  nominal  trajectory,  it  is  assumed  that  the  aircraft 
is  in  a  steady  flight  with  a  constant  velocity  Vj  along  a  straight  path  with  a 
flight-path  angle  y,  and  altitude  hQ  >  2hpUwhere  h„u  is  the  pull-up  altitude  as 
shown  in  Figure  34.  This  type  of  flight  is  assumed  to  be  accomplished  by 
adjusting  the  elevator  angle  (i.  e. ,  stabilator  in  F-4  case)  and  the  magnitude 
of  the  engine  thrust. 


Figure  34.  Nominal  Dive  and  Pull-Up  Trajectory 


In  thjg  second  part,  it  is  assumed  that  the  aircraft  flies  with  a  constant 
velocity  V2  and  a  constant  pitch  rate  Q0.  Again  this  type  of  flight  is  assumed 
to  be  accomplished  by  adjusting  the  elevator  angle  and  the  magnitude  of  the 
engine  thrust. 


In  the  following,  the  trim  equations  for  the  dive  phase  of  the  maneuver 
are  developed  first.  Then  the  constant- g  pull-up  case  is  treated. 

Development  of  the  Missing  Nominal  Values  for  a  Dive  Maneuver  —  It  follows 
6rom  the  abcve  descriptions  that  the  mathematical  implications  associated 
with  the  given  conditions  are: 


Constant  Speed  =»  V  = 


(7.8) 


Steady  Flight  =»  I  v  =  0  and  q 


(7.9) 


Symmetric  Flight 
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(7. 10) 


Z  +  ZT  +  mg  cos  0  =  mw  (7.16) 

L  =  p  s  0  (7.17) 

M  +  Mt  *  Iyq  (7.18) 

N  s  r  =  0  (7.19) 
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0=0  or  equivalently  /  1  =  0 


(7.20) 


Therefore  the  problem  is  to  find  0(t),  6  (t),  and  T(t)  such  that 


V  -^np 

_+_  -gsine  =  0 


(7.21) 


^  +t=?  +  g  cos  0  =  0 
m  m  ° 


(7.22) 


M  ,  MT 
I  I 

y  y 


(7.23) 


Theoretically  speaking,  it  is  impossible  to  produce  a  steady  flight  during  a 
dive  motion  due  to  changes  in  the  altitude  and  Mach  number  parameters. 
However,  for  all  practical  purposes  the  contributions  of  u,  w  and  q  terms 
are  negligible. 

It  should  be  noted  here  that  in  the  solutions  of  the  above  equations,  an 
equivalent  state  a  is  used  instead  of  0,  since  the  aerodynamic  forces  and 
moments  are  functions  of  this  variable.  After  having  found  the  value  of  a, 
0(t)  and  T(t)  are  computed  from 


0  =  a  +y 


(7.24) 


cos  0/2 


sin  0/2 


(7.25) 
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i&m&K's&r'nsr-? 


N  =  r  =  0 
0  =  q 
0  =  0 
i  =  0 

xg  =  ( cos  6)  u  +(sin  6)  w 


ye  -  o 


-(sin  0)u  +(cos  $w 


(7.35) 
(7.  36) 
(7.  37) 
(7.  38) 
(7.39) 


(7.40) 


(7.41) 


Now,  substituting  u  =  V  cos  a,  w  =  V  sin  a  into  (7.30)  and  (7.  32)  yields 
X 

+  — _  «  Sin  q  _  qv  sin  a  =  u  < 

mm 

Z 

~  +  — i-  +  g  cos  0  +  qV  cos  a  =  w  ' 

mm 


(7.42) 


(7.43) 


-  4 

y  y 


(7.44) 


which  shows  that  the  set  (7.42),  (7.43',  and  (7.44)  reduces  to  set  (7.21), 
(7.22),  and  (7.23)  when  q  is  set  to  zero.  Here  u  and  v  are  functions  of  t 
and  u  f  0,  w  /  0.  If  contributions  due  to  u  and  w  are  neglected  (this  is  a 
common  procedure  in  algebraic  trim),  the  set  of  equations  to  be  solved  be¬ 
comes  almost  the  same  for  the  two  parts  of  the  nominal  trajectory. 

In  practice,  usually,  pull-up  maneuver  is  defined  by  specifying  the  so- 
called  load  factor.  Assuming  w  **  0  in  (7. 32 ),  one  can  write 


q  =  JL 
u 


The  load  factor  is  defined  as 


-(Z  +  ZT) 


-  coo  9 


(7.45) 


n  = 
z 


“(Z  +  zT) 


(7.46) 
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In  term.3  of  this  quantity  the  q  terms  in  (7.42)  and  (7.43)  become 


qV  sin  a  =  g(n  -  cos  6)  tan  a 
z 

(7.  47a) 

qV  cos  a  -  g(n  -  cos  6) 
z 

(7.  47b) 

Now,  for  a  symmetric  flight,  the  engine  thrust  variables  are  given  by 

II 

T(cos  ?  cos  ri) 

(7.48) 

ZT 

-  T(cos  ?  sin  r\) 

(7.49) 

and 

mt  = 

=  X,jAz“  Z^Ax  =  T(Az  cos  %  cos  ti  +  Ax  cos  ?  sin  ri) 

(7.  50) 

Observation  of  (7.48)-(7.  50)  togetlier  with  (7.  21)- (7. 23),  or,  equivalently, 
(7.42)-(7.44  ,  indicates  that  the  total  thrust,  T,  enters  into  the  equations 
linearly.  That  is. 

• 

u  = 

T  X 

ci  g  sin  0  -  qV  sin  a  =  0 

lmm  ° 

(7.  51) 

W  = 

T  Z 

c„  —  +  —  +  g  cos  e  +  qV  cos  a  =  0 
<2  m  m 

(7.  52) 

and 

• 

q  = 

o 

CO 

vTM 

ii 

o 

(7.  53) 

where 

II 

O 

COS  §  COS  T| 

(7.  54) 

C2  = 

-  cs?  sin  ri 

(7.  55) 

ca  = 

(Az  cos  §  cos  r\  +  Ax  cos  %  sin  ti) 

(7.  56) 

For  this  reason  the  number  of  equations  to  be  solved  can  be  reduced  to  two 
by  eliminating  T. 
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-f^r, — ;■.  *-«!'^f7 


Solving  (7.  51)  for  T  yields 

T  =  —  [mg  sin  +  mqV  sin  a  -  X] 

C1 

and  substituting  this  into  (7.  52)  and  (7.  53)  one  obtains 
w  =  g(cos  0  -  c4  sin  0)  +  qV(cos  a  -  sin  a) 

+  C(c4X  +  Z  >/m]  =  0 

q  =  [c  {mg  sin  0  +  mqV  sin  a  -  X)  +  M]  /I  =  0 
o  y 

where 

Cj  =  cos  ?  cos  q 
c^  =  tan  q 
c5  =  (Ay  -  c4Ax) 


{7.  57) 


{7.  58) 

(7.  59) 


(7.  60) 


Hence  the  problem  is  reduced  to  finding  a  and  6S  such  that  (7.  58)  and  (7.  59) 
are  both  zero. 

To  summarize,  in  the  development  of  the  trim  equations  it  is  assumed 
that  the  adjustable  parameters  are  the  magnitude  of  the  engine  thrust  and  the 
elevator  deflection  angle: 

•  Trim  Equations  - 


w  =  g(cos  0  -  c4  sin  0)  +  qV(cos  a  -  sin  a) 

(c  X  +  Z) 

+  — - - 

m 

q  =  [M  -  c.  (  X  -  mg  sin  0  -  mqV  sin  a)]/ 1  =  0 

where,  neglecting  a 

y  =  yQ  +  qAT 

0  =  y  +  a 


(7.58a) 

(7.  59a) 


(7.61) 
(7.  62) 
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and 


M  sm  a 


[n  -  cos  0]  pull-up 


Trim  Value  of  Engine  Thrust 

T  =  (mg  sin  6  +  mqV  sin  a  -  X  )/c^ 
Trimmed  Accelerations 


xT+  X 


-  g  sin  6  -  qV  sin  a 


(7.47a) 


(7.57a) 


(7.63) 


Z  +Z 

w  =  - +  g  cos  6  +  qV  cos  a 

m  6 


q  =  <M  +  MT)/Iyy 


(7.64) 

(7.65) 


In  the  following,  a  procedure  is  developed  by  which  0(t),  6s(t)  and  T(t) 
are  found  at  time  points  [tk] ,  k  =  1, 2, -  The  method  is  based  on  the  con¬ 

cept  of  "finite  number  iterations." 


Development  of  the  Finite  Iteration  Algorithm 


In  short  notation,  (7.  58)  and  (7.  59)  are  expressed  with  the  aid  of  (7.  61) 
and  (7. 62)  as 


f.  (Mach,  h,  a,  6n)  =  0 
J.  s 

f2  (Mach,  h,  a,  dg)  =  0 


(7.  66) 
(7.67) 


where  mach  number  and  altitude,  h,  are  parameters  and  a,  6g  are  the  real 
solutions  of  these  nonlinear  algebraic  equations  [42,  43] . 

In  the  fallowing,  we  shall  formally  develop  a  discrete  version  of  the 
method  given  in  [42],  for  solving  (7.  66)  and  (7.  67)  since  the  partial  deriva¬ 
tives  of  fi  and  fr>  cannot  be  evaluated  analytically  as  assumed  in  [4°] .  During 
this  development,  the  parametric  dependence  of  (mach,h)  will  be  sv  "r,essed 
for  the  writing  ease. 


Define 


f  0 1  (®»  5g»  T1»T2)1  rfi<a»  6s)  *  TlVao’6so* 


.02(o'*  5g,  t1^2)J  l/2<«*  5&)  -  T2f2^o’  6s0> 


(7.  68) 


where  aQ,  6S  are  initial  estimates  for  a  and  respectively,  and  Tj,  t2  are 
unknown  parameters.  Clearly  for  the  initial  vAues  Tj=t2^1 

r^i(v6s  •ti*t2)  ° 

=  (7.  69) 


r0l(ao'  6S  ,T1'T2) 
° 

'o' 

\j>2^a>  6s  ,Tl»T2^  * 

.0. 

Moreover,  the  values  a  and  6g  for  which 


i~01(q',  6s*  o.o)]  r° 


02.a-,  6g,  0, 0)J  Lo 


(7.  70) 


are  the  solutions  of  (7.  66)  and  (7.67).  Hence,  one  can  find  the  solution  a,  6S 
by  gradually  decreasing  Tj  and  t2  to  zero  from  the  starting  values  Tj  =  t2  =1, 
and  at  the  same  time  determining  a  and  5g  so  that  0j  =  02  =0  for  every  value 
of  Tj  and  t2. 

Now  let  ATj  and  ATo  be  small  perturbations  about  the  initial  values  t  . 
and  To.  Then  by  expantnng  0j  and  02  about  the  initial  point  (a,  6S» ti,t2)  «  is 
possible  +o  find  small  perturbations  Aar  and  A6S  such  that 

.fill  r01(«  +  A«,  6S+A6S,  Tj  +  ATj,  t2+at2)1  f0~| 


L02J  U>2{q  +  Ao'»  6s  +  A6s’  T1  +  AT1'  T 2  +  AT2)  J  L°J 

Ignoring  the  second-  and  higher-order  terms  and  making  use  of  (7.69)  one 
obtains 


A® +  Sg  A6s  +  fl(ao'6s  >ATi  =  0  (7*71) 


rwB  a.  6 
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‘fic-s  '>p»7  *s'<a79,Tcirj*sf,wcr;T 


*•' 


I 


A 


at, 


ha 


a> 


A  a  + 


at. 


a  6, 


a,  6, 


A5s  +  f2(V  6s, 


)AT.  =0  (7.72) 


where  indicated  partials  are  obtained  by  divided  differences.  Now  if  these 
two  equations  are  independent,  then  A a  and  A6a  can  be  obtained  in  terms  of 
given  At  j  and  Atg. 

By  continuing  this  procedure,  the  parameters  and  T2  can  be  reduced  to 
zero  in  a  finite  number  of  steps  and  the  values  of  a  and  $a  corresponding  to 
Tj  =  T9  =  0  become  the  real  solutions  of  (7.  66)  and  (7.  67).  From  the  above 
formal  treatment,  a  numerical  algorithm  for  finding  the  values  of  a  and  6g 
emerges.  The  subroutine  which  implements  the  procedure  is  called  ** 
SUBROUTINE  NOMK. 

In  the  following,  finite  iteration  algorithm  is  compared  with  the  Newton- 
Raphson  process  with  respect  to  convergence. 


} 


t 

i 

\ 


Finite  Iteration  Algorithm  versus  Newton-Raphson  Process 

The  development  given  above  is  closely  tuned  for  solving  a  specific  trim 
problem.  However,  the  method  of  solution  applies  equally  well  to  solving 
nonlinear  vector  equations  of  the  form 

f(x)  =  0  (7.73) 


with  a  vector  argument.  In  this  case,  one  obtains  the  following  vector  equa¬ 
tion  with  initial  guess  solution  xt°); 


Ax  4-  f(x^°H  At^  -  0 


(7.74) 


Then  for  given  value  of  At 


Ax  is  found  as 


Ax 


bf(x) 

dx 


«x(0))  AT(o) 


(7.75) 


Provided  that  the  inverse  exists 
variables  are  updated  by 

T  =  T  -  At^ 

x  =  x  +  Ax 


Starting  with  x  =  x 


(o) 


and  t 


=  1, 


the 


(7.76) 
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and  a  new  Ax  is  computed  from  (7.  75).  This  process  continues  until  t  =  0, 
and  at  that  point  x'*'  =  x  is  the  new  solution. 

To  improve  x^1),  one  may  repeat  the  computations  with  the  new  initial 
vector  f(x(l))  and  with  possibly  larger  value  of  At(1),  as  llustrated  in  Fig¬ 
ure  35.  Note  that  as  AT(i)  -•  1  this  algorithm  reduces  to  the  Newton-Raphson 
process  (Figure  36)  in  which  the  increments  are  computed  [59]  as. 


«*k> 


(7.77) 


and  solution  is  corrected  as 


(7.78) 


Note  that,  in  the  finite  iteration  algorithm,  the  solution  is  approached  in  small 
steps  (i.  e. ,  At  «  1)  and  corrected  in  small  steps.  Consequently,  it  is  not  so 
sensitive  to  the  convergence  problems  as  the  Newton-Raphson  process. 


TRIMMING  WITH  AN  AUTOPILOT 

Figure  37  shows  a  trimming  process  by  an  autopilot.  During  the  nonlinear 
simulation  of  aircraft,  the  error  signal,  (i.  e. ,  the  difference  between  actual 
and  desired  trajectopr  variable)  is  used  in  a  simple  autopilot  to  generate  a  con¬ 
trol  input.  If  the  gains  are  properly  chosen,  e(t)  can  be  maintained  reason¬ 
ably  small  while  obtaining  trim  profile.  The  controller  equation  used  for 
dive  and  pull-up  is  in  the  form: 

VW  =  W +  V-yo  -  ’,<tkn  +  Xy  <7- 79) 

+  Kq[q°  ‘  q<tk)]  +  Kq  qV 

This  equation  is  implemented  as  subroutine  PILOT  in  ADAPS. 


Figure  36.  Newton-Raphson.  Process 


Figure  37.  Trim  by  Autopilot 


SECTION  vm 

DEVELOPMENT  OF  A  PERFORMANCE  MEASURE  FOR 
OPTIMAL  WEAPON  DELIVERY  CONTROLLER  DESIGN 

In  this  section,  a  methodology  is  developed  for  determining  weighting 
matrices  for  optimal  weapon  delivery  controller  design.  Half  probability 
area,  HPA,  is  chosen  as  a  measure  of  weapon  delivery  performance,  and 
the  effects  of  fi'ght  control  parameters,  airframe  dynamics,  measurement 
errors,  and  gust  disturbances  are  related  to  this  measure  by  using  the  over¬ 
all  system  model. 

HPA  is  the  area  of  a  circle  centered  at  the  mean  impact  point  with  0.  5 
hit  probability.  CEP  is  its  radius. (see  Figure  45).  For  normal  distribu¬ 
tions  with  small  cross  correlations,  this  area  can  be  closely  approximated 
in  terms  of  the  impact  covariance  matrix  of  the  bomb. 

For  this  reason,  performance  analysis  of  a  weapon  delivery  process 
can  be  reduced  to  linear  covariance  analysis.  In  the  following,  a  model  is 
developed  for  determining  the  initial  perturbation  state  of  bomb,  and  propa¬ 
gating  the  release  point  errors  to  impact.  Next  the  expression  for  the  HPA 
performance  measure  is  developed  and  its  approximation  in  terms  of  quad¬ 
ratic  cost  is  given.  The  analysis  is  implemented  as  subprogram  PERK 
which  propagates  the  release  errors  to  impact. 


DEVELOPMENT  OF  THE  INITIAL  STATE  OF  BOMB 

The  statistical  description  of  the  perturbation  state  of  bomb  just  after 
release  is  needed  to  propagate  release-point  errors  to  impact. 

The  state  of  bomb  and  the  state  of  aircraft  which  carries  the  bomb  are 
related  to  each  other  by  a  nonlinear  algebraic  equation 

xb  *  \  |x'  £b'  Arb)  (8.1) 

where 

x  *  state  of  aircraft 

E.  «  bomb  orientation  matrix  (transformation  from 
aircraft  to  bomb  axes) 

Arb  *  bomb  position  vector 
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The  perturbation  state  of  the  bomb  at  t*tr  for  fixed  Eb  and  Ar^ 
is  given  by 


5xbjtr)  e  Hhfix(tr 


(8.2) 


where 


(8.3) 


In  this  work  it  will  be  assumed  +hat 


(8.4) 


(i.  e.,  the  bomb  is  at  the  eg  of  the  aircraft  and  oriented  parallel  to  the  aircraft 
axes).  However,  for  completeness,  the  bomb  station  geometry  will  be 
given  briefly  in  what  follows. 


Nominal  State  of  Bomb  at  Release 

The  twelve  components  of  the  state  of  the  plant  undergo  a  jump  at  the 
release  time  (Figure  38).  This  is  due  to  changing  the  plants,  namely  the 
transition  of  dynamics  from  the  aircraft  to  the  weapon  and  adding  ejection 
velocities. 

x„  «♦  x. 

P  b 

To  compute  the  nominal  state  of  the  bomb  at  release  from  the  nominal 
state  of  the  airplane,  the  bomb  station  geometry  is  considered. 

It  is  assumed  that  the  mass  center  Ob  of  bomb  is  located  by  a  vector 
2krb  from  the  mass  center  0  of  the  aircraft  (Figure  39).  It  is  further 
assumed  that  0^  are  fixed  azimuth  and  elevation  angles  from  aircraft  to 
bomb  axes. 

Then,  as  developed  later  in  this  section,  the  nominal  state  of  the  bomb 
at  release  is  obtained  as  follows  (see  Figure  38): 


Transition  of  the  Linear  Velocity  State  - 
X1(W  *  Eb(6b*  ^  {Es  (tr>  x!s  <W  +  WAr} 


(8.5) 


where 
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Figure  40.  Linearization  About  the  Nominal  Release  Point  &c(t  ),  t  ] 


Equations  (8.17)  and  (8.18)  show  that  the  statistics  of  6tj.  (i.e.,  E { 6tr } , 
E  {fitr^  })  must  be  known  to  compute  release  point  error  statistics. 

The  release  time  error  5tr  consists  of  two  components:  (a)  timing  error 
6tj  due  to  delays  in  the  release  mechanism,  and  (b)  timing  error  5tr  in  com¬ 
puting  the  release  time. 

For  closed-loop  fire  control  algorithms  in  which  release  time  continu¬ 
ously  depend  upon  the  current  state,  the  release  time  perturbation  6tr  is  a 


deterministic  function  of  the  state  at  the  nominal  release  time  tr.  That  is 


tr  =  h[x(t)] 


(8.19) 


For  small  perturbations  in  state.  Equation  (8. 19)  yields 

6t  =  r\  '  5x(t  ) 
r  r 


(8.20) 


where 


.  ilL 


(8.21) 


Equation  (8.20)  is  the  linearized  model  of  the  release  computer.  If  instead 
of  6x,  its  estimate  is  available,  then 


6tr  =  r\'  [6x(tr)  +  e(tr)3 


(8.22) 


where  e(tr)  is  the  error  in  the  state  estimation,  and  6x(tj.)  is  the  optimal 
estimate  of  6x(t). 

Thus,  (8.16)  and  (8.22)  define  the  .  slease-point  error  in  terms  of  the 
perturbations  about  the  release  point. 


Perturbation  State  of  Bomb  at  Release 

The  transition  from  the  aircraft  perturbation  state  to  the  bomb  perturba¬ 
tion  state  is  accomplished  using  equation  (8.  2).  The  nonsingular  matrix  Hb 
is  obtained  from  the  linearization  of  (8. 5)-(8. 8)  at  nominal  release  point 
Cx(tr),  tr].  Substituting  (8.16)  into  (8.2)  yields  the  perturbation  state  of  the 
bomb  at  release: 

6xb(tr)  «  Hb  (6x(tr)  +  f  [x^p,  u(tr)]  6tr3  (8.23) 

As  indicated  before,  in  this  work  H.  =  I  will  be  utilized. 
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Perturbation  State  Transition  During  Release  Transient 

Modeling  of  the  state  transition  during  the  release  transient  phase 
(Figure  41)  is  outside  the  scope  of  this  work.  A  detailed  model  is  currently 
being  developed  in  £9, 1 0  3. 

Since  it  is  difficult  to  predict  the  bomb  aerodynamics  (i.e.,  aerodynamic 
forces  acting  on  bomb)  in  the  region  of  wing,  the  state  transition  during  this 
phase  is  taken  into  account  in  a  c'.rnplified  manner  by  introducing  an  additive 
release  transient  error  ?r  with  known  statistics.  (Ejection  velocity  is  in¬ 
cluded  in  §  . ) 


5W  *  6VV+Hr?r 


(8.24) 


where  H  is  the  release  transient  input  matrix.  Thus,  perturbation  state  of 
bomb  just  after  release  is  given  by 

6xb(tr+)  =  Hk  [  6x(t^)  +  f„6t„  3  +  H„ 


r  r  r  t* 


Figure  42  shows  its  structure. 

Statistical  description  of  perturbation  state  of  bomb  at  release,  (mean 
and  covariance),  are  given  respectively  as 


5x.  =  H.  [fix  +  f  6t  3  +  H  C 

do  D  r  r  r  r  r 


(8.25) 


X.  =  H.Cx  +  f  f'  6tJ  iH'  +  H  S H  7 
bo  brrr  r  b  rrr 


(8.26) 


where  Xr  andZ^,  denote  the  covariance  matrices  of  5xr  and  ?r  respectively. 


In  these  equations  <5xr  and  Xr  are  obtained  by  using  the  linear  equations 
of  the  aircraft  together  with  the  deterministic  and  stochastic  inputs  (i.e., 
mean  wind  and  gust). 

The  estimate  of  Str  and  5tr^  involves  basically,  determining:  (a)  the 
contribution  of  the  release  computer  and  (b)  the  contribution  of  uncorrected 
delays  in  release  mechanism. 

The  estimate  of  5r  and  2r  involves  the  knowledge  of  ejection  velocity, 
its  uncertainties  as  well  as  transient  effects.  In  this  work  tney  are  con¬ 
sidered  to  be  arbitrary  input  parameters. 


! 

f 


Perturbation  State  of  Bomb  at  Impact  Planes 


Consider  the  nominal  and  perturbed  impact  points  represented  by 
Cx(tf),  tf]  and  Cx(?f),'tf]  respectively  (Figure  43). 

Analogous  to  (8. 16)  the  impact  error  can  be  written  as 

6x(tf)  =  5x(t^)  +  f[x(tp,  u(tf)]  6tj  (8.27) 

where  f  is  the  derivative  of  the  nonlinear  bomb  dynamics.  In  (8.27)  6tf  is 
computed  using  horizontal  or  vertical  impact  planes.  Horizontal  impact 
occurs  when  the  end  point  of  the  perturb  trajectory  lies  in  the  xy  plane 
(5h  =  0).  Similarly  vertical  impact  occurs  when  the  end  point  lies  in  the 
yh  plane,  (Sxe  =  0). 

Tnese  take  place  v/hen 


5h(tf) 
Hfh  =  "  fi(tf) 


6x  (tf) 
fit  —  ..  ..  -J5-— JL. 

5tfv  -  ie(tt) 


(8.28) 


Substituting  (8.28)  into  (8.27)  yields  the  horizontal  and  vertical  impact  errors 
in  terms  of  the  nominal  impact  error: 


6  xh(tf)  =  Hh  5x(tf) 
6  xy(tf)  =  Hy  6x(tf) 


(8.29) 


Figure  43.  Linearization  About  the  Nominal  Impact  Point  [x(t-),  tf] 


.  *  ;i  ,1  iv -V«^^V®,Sfiff5^*r£Br?*' "XV  T*zxr  Wrf?*S-3?r^l 


where 


H.  = 


f(tJ^_iL 


chf<V  J 


H 


I  - 


f(tf  )c  x 


r_^r 


'X  ~'  f 
e 


where  c^  and  cx  are  vectors  which  pick  up  h  and  xfi  components  out  of  full 
state  vector  x.  e 


From  (8.29)  the  mean  error  and  the  mean  square  error  matrix  can 
readily  be  obtained 


fixh(tf)  =  Hh  6x(tf) 

Xh<tf)  =  Hh  X(tj)  Hh ' 

(Similar  expressions  apply  to  the  vertical  impact  errors. ) 


(8.31) 


DEVELOPMENT  OF  STATISTICAL  PERFORMANC  E  MEASURE 
FOR  WEAPON  DELIVERY  PROCESSES 


For  a  linear  system  driven  by  a  Gaussian  white  noise,  the  mean  and  the 
covariance  of  the  state  are  described  respectively  by 


x  =  Fx  +  GjU  +  G2w 


y  "  Hx 


(8.32) 

(8.33) 


and 


r>j  />j 


$  =  FX  +  XF  '  +  GjUGj'+G^Gg' 


Y  =  HXH 


(8.34) 

(8.35) 


where 


x  «*  E{x},  u  =  E{u},  w  *  E  lw],  y  *  E  iy ), 

X«=  E{(x-x)(x-x)/},  U  *  E{(u-u)(u-u)/ },  Y  =E  {(y-yKy-y)'}  and 
E{ww'} 


(8.36) 
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Since  the  process  is  linear  the  state  turns  out  to  be  Gaussian  also.  For  this 
reason,  the  mean  and  the  covariance  (i  e.,  first  and  second  moments)  com¬ 
pletely  determine  the  statistical  behavior  of  the  state,  and  the  output.  The 
probability  density  function  of  the  state  vector  x  of  size  nx  1  is  given  [47] 

by 


f(x,  t)  = 


-  i  (x*  x) 7  X_1(t)(x-x) 


(2TT)n  det  X(t) 


(8.36a) 


The  loci  of  constant  probability  densities  at  each  time  instant  are  families 
of  concentric  quadratic  surfaces  centered  at  x  and  are  described  by 


(x-x) 7  X  X(x-x)  =  c 

The  density  is  maximum  at  its  center,  and  it  equals 


(8.36b) 


f(x,  t)  = 


V  (2u)n 


(8.36c) 


det  X(t) 


For  weapon  delivery  performance,  the  distribution  of  X4  (i.e.,  position 
coordinates)  are  needed  for  computing  the  probability  of  an  event  described  by 
1x4  eDJ.  In  general  the  probability  density  function  of  X4  in  subspace 
R3,  (i.e.,  marginal  density)  is  obtained  by  integrating  the  joint  density  fx : 


+  00  +eo  +w 

f  (x4>  =  /  /  /  f  (xi,x9,x0,x4)dx  dx  dx 

x4  *  J  J  J  x  *  &  *  4  123 


(1.36d) 


_  00  aco 


However,  for  the  case  treated  here  this  integral  can  be  eliminated  since  the 
output  is  linear  function  of  state  and  is  therefore,  normal  (i.e.,  Gaussian). 
It  follows  that  X4  has  a  mean  and  covariance  given  respectively  by 


x4  =  Hx 


X.  =  HXH 
4 


(3.37a) 


(8.37b) 


where  H  is  a  matrix  of  size  3  x  n  which  selects  the  position  coordinates 
from  the  full  state  vector  x  of  size  nxl.  Thus  its  density  is  described  as 


fx  <x4>  =  r=— 

*  V  (2  TT  V  det  X4 


rk  (x4-i4) '  X4’i(x4-i4) 


u*^S6Sfiai%£ifi Im  i' a.-,*,--  -  -  - 


(8.38) 


This  shows  that  when  components  Xj  are  jointly  normal,  they  are  also 
marginally  normal. 

3 

Now,  consider  a  region  O  in  R  ,  and  an  event  described  by 

{x4cd}  (8.39) 


The  probability  of  occurrence  of  this  event  is 
P(x4ED)  =  /  fx  (x^,  t)  dx  dy  dh 


(8.40) 


where  the  integration  extends  over  the  region  D. 

Region  D  may  be  a  cube  or  a  ball.  For  weapon  delivery  performance 
against  a  specific  target  of  dimensions  a  x  b  x  c,  hit  probability  Pjj  given 
by  (8.40)  is  used  with 


a  ^  ,  a  b  <•  b  c  c 

=  -5  -o  sy 


(8.41) 


If  the  shape  of  target  is  not  specified,  n-dimensional  ball  (n  =  1,  2,  3)  centered 
at  the  mean  ‘  xsed  most  often  for  region  D.  The  probability  that  the  weapon  is 
within  this  bs  time  t  is  also  given  by  (8.40)  with  integration  extending 
over  the  ball: 

D:  (x-i)2  +  (y-y)2  +.(h-h)2  *  (R)2  (8.41a) 

The  radius  R  for  which  P  =  is  referred  to  as  the  "spherical  probable 
error"  and  is  denoted  by  SEP  or  SPE. 


When  the  region  D  is  circular,  the  radius  of  its  boundary  is  referred  to 
as  CEP  or  CPE.  (These  are  confusing  names  attached  to  a  radius.) 

Figure  44  demonstrates  the  evolution  of  the  spherical  region  D  as  a 
function  of  time. 


In  air-to-air  weapon  delivery,  SEP  can  be  used  to  measure  the  delivery 
performance.  For  air-to-ground  delivery,  a  simpler  measure,  CEP  can  be 
used  effectively. 


For  horizontal  targets  CEPjj 


f  v(x,y)  dx  dy  - 


1 

2 


is  obtained  from 


(8.42) 
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Figure  45.  Regions  for  Horizontal  and  Vertical  CEP 
Evaluations 


and 

J  fh(h)  dh  =  J,  Dh:  |h-h|  s  EEP  (8.48) 

DH 

Table  III  shows  the  terminology  associated  with  the  radius  of  the  half 
probability  ball  for  one-,  two-  and  three-dimensional  balls. 


Table  III.  Terminology  for  the  Radius  of  Half  Probability  Ball 


Dimension 
of  Ball 
(n) 

Symbol  for  the 
Radius  of  Half 
Probability  Ball 

Name  for  the 

Radius  of  Half 
Probability  Ball 

3 

SEP 

Spherical  Error  Probable 

2 

CEP,  (CEPh,  CEPy) 

Circular  Error  Probable 
(horizontal,  vertical) 

1 

LEP,  (REP,  DEP,  EEP) 

Linear  Error  Probable 
(range,  deflection,  elevation) 

In  general,  the  integrals  given  by  (8.40),  (8.41a)  and  (8.43)  cannot  be 
evaluated  analytically.  The  results  from  a  computer  solution  for  the  case 
where  x  and  y  are  independent  stochastic  cesses  is  given  in  [44]. 
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In  [45,  46],  CEP  is  expressed  as 


CEP  =  Q(p)  a  .  (8.49) 

and  this  dependence  is  exhibited  as  shown  in  Figure  46,  where  Q(p)  is  referred 
to  as  the  CEP  parameter  and  is  a  function  of  the  ratio 

G. 

OSp  =  1  (8.  50) 

and  a.  greater  of  (ax»ay}»  ^  =  smaller  of  {ax>  c^},  and  ox,  are  the  cross¬ 
range  and  down-range  variances,  respectively. 


i 

3 

a 


Figure  46.  CEP  Parameter  versus  Standard  Deviation  Ratio 


Now,  we  shall  introduce  a  modified  performance  measure,  "Half- Proba¬ 
bility-Area  (HPA)"  defined  as 

HPA  =  (n)  (CEP)2  (8.  Jl) 

and  relate  this  to  the  quadratic  measure  described  by 

J  =  tr  (HX(t)H  0  (8.  52) 
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where  tr  is  the  trace  operator  and  H  is  the  output  matrix  which  selects  and 
weighs  appropriate  elements  (i.  e„  a^,  cR,  a£)  of  the  covariance  matrix.  To 
be  specific,  we  wish  to  find  weighting  parameters  a  and  0  such  that  the 
quadratic  cost  given  by 

J  =  <tt>  (a)  [1+  (p$20]  a2  (8.53) 

where 

a-i 

p  =  J-i  1  (8.  54) 

l 

is  a  good  approximation  to  HPA. 

For  a  close  fit,  a  and  0  should  be  functions  of  p.  For  simplicity,  we  shall 
assume  constant  a  and  0  in  the  interval 

0  £  p  £  1. 

Approximation  to  Weapon  D~ livery  Performance 

The  CEP  parameter  Q(pJ  is  approximated  by  two  line  segments 

Q(p)  —  k.  (1  +  p-p)  i  =  1,  2  (8.55) 

for  the  ranges  0  s  p  5  0.  2  and  0.  2  £  p  £  1.  0,  as  shown  by  dotted  lines  in  Fig¬ 
ure  46.  The  values  for  k.  and  p^  are: 

0  £  P  £0.2  0.2  £  p£  1.0 

kx  =  0.  6744  k2  =  0.  585 

Pj  =  0.  24  P2  =  1 

Figure  47  shows  the  HPA  parameter  §(p)  defined  by 
Q(p)  =  [Q(p)3 2 

A 

Also  shown  is  an  almost-bound  to  Q(p)  defined  by 
Q(p)  =  or [1  +  0  p2] 

where  a  =  j  and  0  =  2. 

With  these  values,  the  quadratic  cost  becomes 


(8.  56) 

(8.57) 

(8.58) 


141 


mrm 

ihhbs 

HB 


EWIV^. 


^nfr 


2 1  2 


J  =  -g-  [1  +  2p  ] 


(8.  59) 


j  •  «i °?  +  <52-  °i><5 


<8.  60) 


where 


q.  -  -p  and  q  •  ~  tt 


(8.  61) 


Finally  from  this  quadratic  cost  (i.  e.  half- probability  area)  approxi¬ 
mate  CEP  can  be  computed  as 


CEP 


(8.62) 


Development  of  Design  Performance  Index 

It  was  shown  in  the  previous  section  that  for  normal  distributions  with 
small  cross-correlations,  the  half- probability  area  HPA  is  almost  bounded 
by  the  quadratic  cost: 


J  *  0i2  +  qj  oj2 


(8.63) 


where 


cjx  =  downrange  standard  deviation 
ct  =  crossrange  standard  deviation 

y 

or.  =  greater  of  £ax, 
a  =  smaller  of  (ox,  cy] 
q.  =  rr/2 

qj  =  " 

Equation  (8. 63)  can  be  written  in  terms  of  the  covariance  matrix  as 


HPA  a  tr  [HX(t)H')  =  tr  £X(t)QQ} 


(8.  64) 


where 


which  selects  and  weighs  the  downrange  and  cross -range  variances  out  of  the 
full  weapon  state  covariance  matrix. 


DEVELOPMENT  OF  THE  RELEASE  ERROR 
PROPAGATION  PROCESS 

In  this  subsection  the  bomb  covariance,  CEP  measures  and  the  equivalent 
quadratic  weighting  matrices  are  developed  in  terms  of  the  release  state 
covariance,  using  the  results  of  the  initial-state  development  and  the  wind 
model.  This  analysis  is  implemented  in  subprogram  PERK. 

Figure  48  shows  the  vector  state  diagram  of  a  linearized  bomb  model 
with  a  wind  driver. 


Figure  48.  Bomb  Dynamics  with  a  Wind  Driver 


The  state  equations  are  given  by 

x  =  F  x  +  G  n  (8. 66) 

w  w  w  w  Mw  ' 

ro 

wXw  (8.67) 
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(8.  68) 


— 'T^p 3T^j4>rv^  ^%’ 


and 


x. 


=  Fbxb  + 


Gblub  + 


Gb2^w  +  Gb3^ 


w 


n> =  Hbxb 

letting  x  ■  col  (xb|xw)  and  substituting  (8.67)  into  (8.68)  yields 
x  =  Px+GjUh+G^+Qjt, 
yb=  Hx 


(8.69) 


(8.  70) 
(8.71) 


where 


F  = 

'Fb 

(Gb3>(Hw>l 

G  -|Gbl)  G  - 

1  Gfc>2 1  r  #  0  \ 

X 

6 

F 

w 

’  Gi  "inn  ’  g2  - 

1  0  !■  °3-lGw| 

and  H  =  (Hb  |  0) 


(8.  72) 


From  (8.  70),  (8.  71)  and  (8.  72)  the  mean  and  covariance  responses  are 
determined  using  (8.  32)  through  (8.  35). 

For  convenience,  the  total  state  covariance  matrix  is  separated  into  two 
components,  (a)  homogeneous  and  (b)  forced.  This  enables  one  to  evaluate  an 
impact  covariance  matrix  for  arbitrary  initial  covariance  matrices.  Obv. .  «’^ly, 
if  only  one  initial  covariance  matrix  is  considered  it  is  better  to  integrate 
(8.  34)  with  nonzero  initial  conditions  to  obtain  total  covariance. 


It  should  be  noted  that  forced  covariance  component  is  obtained  by  inte¬ 
grating  the  differential  equation  (8.  34)  rather  than  evaluating  the  integral 


0  (t,  s)  W(s)  0 * (t,  s)  (is 


(8.  73) 


Where  0(t,  s)  is  the  state  transition  matrix. 

The  total  covariance  at  any  time  instant  is  given  by 

X(t)  =  0(t,  tr)  Xr+  <f>  (t,  tr)  +  Xf(t),  t  £  tr  (8.  74) 


145 


■  ■■  - . -  mm  - 


nwiliitflnn  --  .  nrliiltiftiiioa 


As  was  shown  previously,  the  quadratic  approximation  to  HPA(t)  is  given 

by 

HPA(t)  2:  J(t)  =  tr  £X<t)Q0)  (8.75) 

Substituting  (8.  74)  into  (8.  75)  yields 
J(t)  =  tr  {(0Xr+0/  +  Xf(t))  Qq} 
or 

J(t)  =  tr  ((X^)  W  Qo0)  +  (Xf(t))  (Qo)}  (8.  76) 

The  matrix  defined  by 

Q  £  0'Qo0  (8.77) 

where  0  *  0(t,  tp)  will  be  referred  10  as  the  ctate  covariance  weighting  matrix 
or  the  propagation  matrix. 

The  approximate  radius  6Sp(t)  is  computed  from 

CEP(t)  t£tf  <8-78> 

As  developed  previously,  the  non- zero  elements  of  Q0;  (qx,  qy)  or 
(qv,  q^)  are  determined  from  (c^2,  ay2)  or  (ay  ,  o\i)  by  a  simple  test. 
Namely,  in  each  pair,  greater  variance  is  associated  with  weighting  value 
it  /  2  and  smaller  with  tt. 

At  impact  it  follows  from  (8.  74)  that 

X(tf)  =  0<tf,tr)Xr+  0#(tftr)  +  Xf(tf)  (8.79) 

Substituting  this  into  (8,  31)  yields  the  impact  covariance 

Xh(tf>  =  Hh  U  #'  +  Xf]  Hh'  (8.  80) 

where  0  and  Xf  are  matrices  evaluated  at  impact  time  tf . 

HPA  at  impact  is  then  given  by 

HPA(tf)  2:  j(tf)  =  tr  {X(tf)Q03  (8.81) 

Substituting  (8. 80)  into  (8.  81 )  yields 
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J(tf)  -  tr  {Hh  (0  X^  0'  +  Xf)  Hh/Q<) } 

(8.  82) 

or 

J(tf)=  tr  ^  C0,Qo0]+XfQo} 

(8.83) 

where 

Qo  ■  V  Qo  Hh 

(8.  84) 

is  called  the  impact  propagation  matrix". 


DEVELOPMENT  OF  THE  VARIANCE  CONTRIBUTION  MATRIX 

The  variance  contribution  matrix,  V,  displays  the  effects  of  the  pertur¬ 
bation  state  vector  components  onto  the  position  error  variances  at  impact. 
Its  brief  development  is  presented  in  what  follows: 

The  diagonal  elements  of  X(tf),  can  be  expressed  as: 


x(tf)  -  ^1*1^  +  ^/2x2r+  +  *  *  •  +  ^/nxnr+  +  xf*V 


(8.85) 


wnere 


x-^  m«>i  xn  are  the  column  vectors  of  the  X  .  matrix,  that  is 
r+  r+ 


"r+  =  Xl^i.  ^X2  .  I  *  *  *  I  **1  1 

r+  r+  r+ J 


(8.  86) 


and  x(tf),  Xf(tj-)  are  vectors  made  up  from  the  diagonal  elements  of  X(tf)  and 
Xf(tf)  matrices,  respectively.  The  elements  of  the  weighting  matrices  and 
covariance  vectors  xj^  for  j  =  1, . . .  n  are  given  by 


Vll  0lj012 

%021  02j022 


•30ln] 


02j02n 


6x^6x^ 


5x26xj 


V  * 


,  x,  =  E 


(8.87) 


Jr+  ?  • 


i 

|^0nj0nl  0nj^i2  .  0nj0r 


6x.6x. 
3  3 


6x„  6x . 


i  . . nw  l  .ww  yppiii 1  .i " 11 1 » 


In  the  decomposition  represented  by  (8.  85),  each  vector  in  the  sum,  e.  g. , 


Vw 


is  identified  as  the  contribution  of  the  perturbation  state  component  6xj  to  the 
impact  covariance  vector.  Note  that  in  this  identification  the  cross-variance 
terms  are  also  included  as  given  by  (8.  87). 


If  position  error  covariance^  (cxe^»  aye^>  0h^)  about  the  nominal  impact 
point  are  of  interest  only,  the  three  rows  of  (8.  85)  need  be  evaluated. 


The  matrix  which  is  made  up  of  the  vectors  t//j  xj^  is  called  the  variance 
contribution  matrix: 


V  = 


1 


=  I  •  •  •  fof 1  KXn^ 


Each  row  sum  of  this  matrix  gives  a  component  of  the  error  variance 
about  the  nominal  impact  point.  Dividing  the  elements  in  each  row  by  their 
row  sum  gives  a  normalized  variance  contribution  matrix,  in  which  each  ele¬ 
ment  shows  the  relative  contributions  of  state  components  onto  position  error 
variances  about  the  impact. 


The  normalized  variance  contribution  matrix  is  used  to  identify  those  state 
components  which  are  important  contributors  to  impact  error  variances. 
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SECTION  IX 


DEVELOPMENT  OF  A  METHOD  FOR 
NONSTATIONARY  OPTIMAL  WEAPON  DELIVERY  CONTROLLER  DESIGN 


In  this  section,  a  method  is  presented  for  designing  an  optimal  non¬ 
stationary  perturbation  coniroller  for  the  precision  weapon  delivery  processes 
[47,  11]. 

The  optimization  problem  considered  here  involves  direct  minimization 
of  the  CEP.  Strictly  speaking,  CEP  is  the  radius  of  a  0.  5  probability  circle 
centered  at  the  mean  impact  point  of  a  bomb  a*  target.  For  normal  distributions 
with  small  cross  correlations,  Hie  area  of  this  circle  can  be  closely 
approximated  in  terms  of  the  impact  covariance  matrix. 

As  shown  in  Section  VIII,  about  a  nominal  trajectory,  the  impact  covariance 
matrix  can  be  expressed  in  terms  of  the  state  covariance  matrix  at  bomb 
release.  By  this  process,  controller  optimization  for  the  precision  bomb 
delivery  is  reduced  to  optimization  with  a  terminal  time  performance  index, 
in  which  the  state  deviations  at  the  nominal  release  time  are  penalized  by  a 
weighting  matrix  which  depends  upon  the  nominal  bomb  trajectory 
sensitivities. 

In  the  following,  the  statement  and  solution  of  the  optimization  problem 
corresponding  to  continuous  time-varying  processes  are  given  first  for 
completeness.  Then  the  discretized  model  and  its  solution  are  developed. 

This  solution  is  implemented  as  a  program  called  DISCOP. 


STATEMENT  OF  THE  PROBLEM 
Given  the  linear  system 

x(t)  =  F(t)x(t)  +  G1(t)u(t)  +  Ggttfv^t)  +  GgWri^t) 

r(t)  =  +  DjftMt)  +  D2(t)vjt)  (9. 1) 

m(t)  =  H2(t)x(t)  +  T)  2(t) 

wher  ■  a  deterministic  (kiiown)  time  function,  the  inputs  T)j(t)  and 
Tl 2^t)  c  .pendent  white  noise  processes 

E  jTljWrijCij)'}  =  WjUWt-tp 

E  {1l2(t)il2(t1>#}  =  W2(t)6(t-t1)  (9.2) 

E  iri1(t>Tl2<t1)'f  =  WgdWt-tj) 
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and  the  initial  state  mean  and  covariance  are  known 


E  |x(o}}  =  x(o)  (9.  3) 

E  |(x(o)  -  x(o»  <x(o)  -  x(o»'}  =  X(o) 

Find  the  linear  control  functional  of  past  and  present  measured  outputs  m(t) 
u(t)  =  L(t,  m[o,  t]),  (9.  4) 

that  minimises  the  quadratic  J 

T 

J  =  tr  rQ(T)S(T)  +  V(T)F'T)  +  J  (Q(t)S(t)  +  V(t)R(t))dtJ  (9.  5) 

o 


where  R(t)  is  the  mean  response  matrix  defined  by 

R(t)  =  r(t)r{t)y  (9.6) 

S(t)  is  the  response  covariance  matrix  defined  by 

S(t)  =  E  {(r(t)  -  r(t))  (r(t)  -  r(t))'}  (9.7) 

and  Q(t),  V(t)  are  the  weighting  matrices,  assumed  to  be  symmetric  and  non 
negative  definite  for  all  te  [0,T],  the  matrices  D/1Q(t)D1  and  D/1V(t)D1  are 
assumed  positive  definite  for  all  te  [0,T] 

Q(T)D1(T)  =  V(T)D1(T)  *  0  (9.8) 


SOLUTION  OF  THE  PROBLEM 

In  this  subsection,  solution  to  the  optimization  problem  posed  above  is 
given  in  terms  of  differential  equation  formulation.  It  is  shown  [46]  (hat  the 
minimization  of  J  can  be  separated  into  a  deterministic  problem  and  a 
stochastic  problem,  and  that  tne  solutions  of  these  two  problems  can  be 
combined  to  form  the  optimum  controller. 


Minimization  of  the  Continuous  Model 

The  J_  minimization  problem  can  be  divided  into  the  control  of  the  mean 
response  r(t)  and  the  control  of  the  response  deviation  from  the  mean, 
r(t)  -  r(t). 
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By  defining  x,  m,  u  to  be  mean  responses,  and 

1  ' 

tf  - 

a  1 

X*  =  X  -  X 

•a  ■- 

m*  =  m - m 

1  \ 

1  ^ 

u*  =  u  -  u 

(9.  9) 

r#  =  r  -  r. 

1  i 

$  : 

■A 

to  be  deviations  from  the  mean,  there  results  the  two  sets  of  system  equations 

i  < 

: 

* 

JS  r- 

i(t)  =  F(tWt)  +  G1<t)u(t)  +  GgW^t) 

U  ) 

1  J 
p  ’ 

r(t)  =  Hx(t)x(t)  +  DjWutt)  f  D^tjv^t) 

(9.  10) 

«  ' 

$  : 

\  < 

m(t)  =  H1(t)x(t) 

1  , 

*  j 

and 

^  ’‘J 

I  1 

x*(t)  =  Flt^tO  +  G^tJu^t)  ^-G^t^t) 

1  j 

1  ' 

i*(t)  =  H^Ox^t)  +  Dx(t)u*(t) 

(9. 11) 

s’  i 

i 

J 

m*(t)  =  H2(t)x*(t)  + T!2(t) 

!  i 
*  : 

The  cost  J  becomes 

N 

i  : 
i  ; 

T 

J=  [r(T)'v(T}T(T)  +  J7(t)/V(t)r(t)dt] 

1  i 
;  -* 

j  i 

0 

T 

(9. 12) 

»  "! 

i 

+  E[r*(T)/Q(T)r*(T)  +  J  r*(t)'Q(t)r*(t)dt] 

o 


Since  u*  in  (9. 11)  does  not  in  any  way  affect  r  in  (9. 10),  and  u  in  (9. 10)  does 
not  in  any  way  affect  r*  in  (9. 11),  the  controls  u  and  u*  rnay  be  designed 
separately  to  minimize  their  respective  contributions  to  J. 


Development  of  the  Optimal  Mean  Control  Law  —  The  mean  control  u  will  be 
determined  first  L47J.  Given  (9. 10),  the  functional  J  where 

T 


J  =  r<T)' V(T)r<T)  +  Jr(t)' V<t)r(t)dt 
is  to  be  minimized. 


(9. 13) 


o 
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Examination  of  the  responses  whose  terminal  behavior  is  to  be  controlled 
reveals  that  the  equations  for  those  responses  do  not  contain  the  final  control 
u(T).  It  can  therefore  be  assumed  that  the  contribution  of  IT(T)  toYCT)'  V(T)r{T) 
is  zero.  That  is 

V(T)Dj\T)  =  0  (9.14) 

The  above  problem  is  then  a  Bolza  variational  problem.  By  writing 

T 

J  =  GjX(T)  +  J  G2(x(t),  u(t),  t)  dt  (9. 15) 


The  Hamiltonian  for  the  problem  is 

H  =  G2  +  X'x  (9. 16) 

and  the  control  u(t)  is  defined  by  the  equations 

*H  ,  0 
du(t) 

4^  =  -  ^(0  (9. 17) 

dx(t) 

dG1 

=  X'(T) 

ax(T) 

With 

hT  =  CHx(t)x(t)  +  Dj(t)u(t)  +  D2(t)vw(t)J'  V(t)[H  X(t)x(t)  +  Dx(t)u(t) 

—  T  t  r-  _  _  _  (9.  18) 

+  D2(t)Vuj(t)J  +  X(t)  CF(t)x(t)  +  G1,.t)u(t)  +  G3<t)v Jt>] 

then 

=  -  X'(t)  =  F(t)/  X  (t)  +  2Hx(t)'  V(t)[H1(t)x(t)  +  Dx(t)u(t) 


=  X'(T)  =  2H1(T)/  V(T)[  HjCT)'  x(T)  +  D2(T)V(j(;(T)] 

=  0  =  G^t)'  X(t)  +  2D1(t)'  V(t)[Hx(t)x(t)  +  Dx(t)u(t)  +  D^Ov^t)] 


iiL 

^x(t) 


dG 


1 


dx(T) 

d  H 
bu(t) 
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By  making  the  additional  substitutions  [49] 

X(t)  =  2[Pv(tRt)  +  g(t)] 

v  (g^  jg) 

u(t)  =  Ky(t)x(t)  +  fy(t), 

and  assuming  that  the  inverse  [Did)7  V(t)Di(t)3“*  exists  for  t  e  [o,  t],  there 
results  the  familiar  Riccati  end  conditions: 

PV(T)  =  H,(T)'  V(T)H.(T) 

V  1  1  (9. 20) 

g(T)  =  Hi(T)/  V(T)D2(T)Vuj(T), 

the  backwards  differential  equations 


-  Py(t)  =  F(t)'  Py(t)  +  Py(t)F(t)  +  H^t)'  V(t)Hj(t) 

-  CP^tjGjd)  +  Hjd)'  V(t)D j(t)3 CD^t)7  VttjDjd)]"1 
.  fc^t)'  Py(t)  +  Dj(t/  V(t)H1(t)  ]  (9. 

-  g(t)  =  F(t)#  g(t)  +  [Pv(t)G2(t)  +  H1(t)/  V(t)D2(t)3vw(t) 

-  [PyttjGjU)  +  HjU)'  V(t)D1(t)][D1(t)/  VCOD^t)]"1 
.  [G^t)'  g(c)  +  D^t)'  VdJD^tVv^d)], 

and  the  controller  equations 

I<v(t)  =  -  [Djft)7  V(t)D1(t)]_1[G1(t)/  Py(t)+  Dx(t)7  VttjH^t)] 
fy(t)  =  -  [D^t/  V(t)D1(t)]"1[Gi(t)7  g(t)  +  D^t)'  V(t)D2(t)vJt)3 

These  three  sets  of  equations  completely  define  the  mean  control  u(t). 

Development  of  the  Optimal  Stochastic  Control  Law  —  Given  the  system 
equations,  (9.  ll),  it  remains  to  find  the  controller 

u#(t)  =  L(t,  m*[o,  t3) 

minimizing  the  quadratic  J* 


(9.21) 


(9.  22) 


J*  =  E[r*(T)7  Q(T)r*(T)  +  J  r*(t)7  Q(t)r*(t)dt3 


(9.23) 


Let  x*(t)  be  the  sum 
x*(t)  =  x^(t)  +  xtt) 


(9.  24) 
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where  x^(t)  is  defined  by  the  orthogonality  condition 

E[x1(t)/  x(t)  jm*(o,  t),  u*(o,  t)]  =  0  (9. 25) 

This  is  the  expected  value  of  the  product  xi(t)'  x(t),  given  present  and  past 
output  measurements  m*(t)  and  past  inputs  u*(t).  Let  Kq(0  be  the  set  of 
gains  defined  by  (9. 20),  (9. 21)  and  (9, 22)  when  V(t)  is  replaced  by  Q(t).  It 
is  asserted  that  the  controller  L(t,  m#(o,  t))  minimizing  die  quadratic  J*  is 

u*(t)  =  KQ(t)Xl(t) 

That  is,  the  optimum  control  input  u*(t)  is  the  product  of  the  state  estimate 
xj(t)  defined  by  the  orthogonality  relation  (9.  25),  and  the  gains  Kq(0  defined 
by  the  Riccati  equations,  (9,20),  (9.21)  and  (9.22). 

This  assertion  is  known  at  the  "separability  property".  It  permits 
separating  the  J*  minimization  into  two  problems: 

•  The  determination  of  the  state  estimate  x*(t) 

•  The  determination  of  the  controller  gains  KQ(t)  that  would  be 
employed  if  the  entire  state  x*(t)  could  be  measured,  and  the 
system  inputs  ^(t)  were  known 

State  Estimation 

The  state  estimate  xj(t)  must  be  generated  to  complete  the  controller 
design.  The  problem  of  generating  state  estimates  xi(t)  satisfying  the 
orthogonality  (9. 25)  has  been  completely  resolved  [50  J.  It  is  well  known 
that  the  orthogonality  condition  (9. 25)  implies  that,  with  x-fc(t)  Gaussian, 
xi(t)  is  the  conditional  expectation 

x^(t)  =  E  |x*(t)  |m*[o,  t],  u*[o,  t]}  (9.26) 

It  can  be  shown  that  x-(t)  could  be  generated  by  a  linear  transformation  of 
m#(t)  and  u*(t) 

x^t)  =  L*  {t,  m*[o,  t]  ,  u[o,  t]  f  (9.27) 

In  [49],  the  appropriate  linear  transformation  L*  is  developed  for  the  case 
where 

rank[W2(t)]  -  dim[m*(t)] 

w3(t)  =  0 

That  is,  every  nonzero  linear  transformation  of  m*(t)  contains  white  noi3e, 
and  the  inputs  hjOlr^t)  in(9.11)  areunccrrelated.  In  [52],  the  latter  restriction 
(W3  =  0)  is  removed. 


Assuming  the  system  is  of  the  form  (9. 11),  x-(t)  could  be  generated  from 
the  linear  system  1 

xx(t)  =  [F(t)  -  L(t)H2(t)]Xl(t)  +  L(t)m*0'  (9.  28) 

where  L(t)  is  the  solution  of  the  forward  Riccati  equation 

P^t)  =  F(t)PT1(t)  +  P^tjFtt)1  +  G2(t)W1(t)G2(t)/ 

-C(Pr|(t)H2(t)/  +  GjMWgft)  ]  (9.  29) 

(W2(t)"\[W3(t)/G1(t)  +  H2(t)Pn(t)3 
Then  the  estimator  gains  are  given  by 

L(t)  =  [PT|(t)H2(t )'  +  GjttjWgCt)]  W2(t)-1 
where  the  initial  conditions  are 

P^o)  =  covCx*(o)x*(o)']  (9.  30) 

Combining  the  Results 

The  optimum  control  is  the  sum  of  the  deterministic  and  stochastic 
solutions 


u(t)  =  u(t)  +  u*(t) 

=  Ky(t)x(t)  +  fy(t)  +  KQ(t)Xl(t) 

=  KQ(t)[Xl(t)  +  x(t)3  +  [Ky(t)  -  KQ(t)3x(t)  +  fy(t) 
=  KQ(t)x(t)  +  f(t) 

where  f(t)  is  the  deterministic  input  given  by 
f(t)  =  [Ky(t)  -  KQ(t)]x(t)  +  fy(t) 
and  x(t)  is  the  conditional  state  expectation 

x(t)  =  xx(t)  +  x(t)  =  E[x(t)|m(o,  t),  u(o,  t),  v^o,  t)J 


(9.  31) 


(9.  32) 


(9.  33) 


where  x  and  x.  are  generated  by  (9. 10)  and  (9, 28),  respectively.  Figure  49 
shows  the  decomposition  of  the  state  vector. 


DISCRETIZATION  OF  THE  CONTINUOUS  PROCESS  AND  OPTIMIZATION 
OF  THE  DISCRETIZED  SYSTEM 

For  computational  purposes  the  differential  equations  given  above  are 
approximated  by  difference  equations.  Also,  the  performance  integral  is 


X 


Figure  49.  Decomposition  of  the  State  Vector 

approximated  by  a  sum.  This  is  called  the  discretization  of  the  continuous 
optimization  problem.  The  discretized  model  is  derived  with  two  constraints 
in  mind.  Both  of  the  constraints  are  based  on  the  desire  to  achieve 
reasonable  computation  time  for  the  optimization  program  and  at  the  same 
time  maintain  a  sufficiently  accurate  approximation  of  the  differential 
equations,  and  the  performance  integral. 

The  simplest  discretization  is  based  on  the  rectangular  rule.  This 
approximation  would  be  sufficiently  accurate  if  At  is  chosen  sufficiently 
small.  But  the  computation  time  is  inversely  proportional  to  A  t;  to  reduce 
computation  time  it  is  desirable  to  choose  At  to  be  as  large  as  possible. 

In  the  following  discretization  based  on  the  rectangular  rule  is  presented 
first.  Then  a  more  accurate  discretization  with  matrix  exponentials  is 
given.  The  latter  discretization  can  be  used  for  the  high-frequency  dynamics 
of  the  overall  system,  and  former  for  the  low-frequency  dynamics  of  the 
system. 


Discretization  by  Rectangular  Rule 

Choose  a  large,  finite  integer  N  and  let 
T 

At  =  ¥• 

Define 

A(n)  =  I  +  At  F(nAt)  H^n)  =  H^nAt) 

Bj(n)  =  At  G^nAt)  D^n)  =  D^nAt) 

B2(n)  =  At  G2(nAt)  D2(n)  =  D2(nAt) 

B3(n)  =  At  G3(nAt)  H2(n)  =  H2(nAt) 

(n)  =  v^  (nA  t) 


(9.  34) 

W^n)  =  Wx(nAt) 

W2(n)  =  W2(nAt) 

W3(n)  =  W3(nAt)  (9.  35) 
Q(n)  =  Q(nAt) 

V(n)  =  V(nAt). 
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The  above  quadratic  problem  may  then  be  restated:  given  the 

x(n+l)  =  A(n)x(n)  +  Bj(n)u(n)  +  B^njv^n)  +  B^nh^n) 
r(n)  =  H1(n)x(n)  +  D^nMn)  +  D^njv^n) 
m(n)  =  H2(n)x(n)  +  ^(n) 
where 

E  ^Ti1(i)-ni(j)/^  =  < At)** 1  W^ije 
E  ^ri2(i)r|2(3)/|  =  (At)"1  W2(i)6  tj 
E  {t|1(i)ri2(j)/}  =  (At)"1  W3(i)6  .. 

6..  =  1,  ^-Oififj 
E  jx(o)[  =  x(o) 

E  |(x(o)  -  x(o))  (x(o)  -  x(o))'[  =  X(o) 
find  the  linear  functional  (i.  e. ,  piecewise  constant  controller) 
N 

u(n)  =  E  L(n,  i)  m(i) 
i=o 

that  minimizes  the  quadratic  J 

r  (N-l) 


where 


R(n)  =  r(n)r  (n)* 


S(n)  =  E(  (r(n)  -  r(n))  (r  (n)  -  r(n))'} 


(9.  36) 


(9.  37) 


(9.  38) 


(9.  39) 


J  =  tr  t(Q(N)S(N)  +  V(N)R(N)  +  E=  At(Q(n)S(n)  +  V(n)R(n))J.  (9.  40) 


(9.41) 


It  is  assumed  that  Q(n)  and  V(n)  are  symmetric  and  nonnegative  definite  for 
n=0, ...N#  and  Di(n)#Q(n)Dj(n)  and  Dj(n)  V(n)Dj(n)  are  positive  definite  for 
n  <  N,  and 


Q(N)D1(N)  =  V(N)DX(N)  =  0 


(9.  42) 


Discretization  by  a  Matrix  Exponential 


For  a  given  value  of  At  a  more  accurate  approximation  to  equation  (9. 1) 
is  given  by  the  sample-data  form 
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c[(k+l)A  t]  =  eAtF<kAt)  |x(kAt)  +  F_1(kAt)  [i  -  e_AtF<kAt)] 
.  [G1(kAt)u(kAt)  +  G2(kAt)7KkAt)  +  Gg(kAt)v  JkAt)]  } 


<9.  43) 


This  form  is  approximate  in  that  the  various  coefficients  are  not  constant  over 
the  At  intervals,  and  the  control  u(t)  is  continuous  and  not  piecewise  constant. 
The  major  disadvantage  of  equation  (9. 43)  is  that  almost  all  of  the  elements 
of  the  coefficient  matrices  are  nonzero,  whereas  in  equation  (9.  35)  the 
majority  of  the  elements  of  the  coefficient  matrices  are  zero.  Computation 
time  increases  at  least  linearly  with  the  number  of  nonzero  elements  [ill 


SOLUTION  FOR  THE  DISCRETIZED  MODEL  WITH  PIECEWISE  CONSTANT 
CONTROLLER 

The  solution  to  the  a’  >ve  discretized  model  problem  follows  that  pre¬ 
sented  for  the  continuous-model  problem.  The  optimum  control  is  of  the 
form 


(9.  45) 


u(n)  =  Kg(n)x(n)  +  [K  y(n)  -  K^WJx  (n)  +  fv(n)  (9.  44) 

where  x(n)  is  the  a  priori  mean  state 

x(n)  =  E  {x(n)|  (9.45) 

and  x(n)  is  the  conditional  estimate 

x(n)  =  E  |x(n)  |  m(o), . . .  m(n),  u(o), . . .  u(n-l),  vjo), . . .  vjn-l)} 

The  gains  K_.(n)  and  input  fv(n)  are  the  solutions  of  the  backwards  difference 
equations  v  v 

Py(N)  =  H^N)'  VtNjH^N)  (9.  46) 

g(N)  =  H1(N)/  V(N)D2(N)^(N)  (9. 47) 

Kv(n)  =  -  [B.(n)'  Pv<r>+l)B1(n)  +  AtD-(n)'  V^D^n)]"1 
V  r  ,  ,  1  -  (9.48) 

.  [B1(n)/  Py(n+l)A(n)  +  AtD^n)  VfajH^n)] 

fv(n)  =  -  iBAn)'  Pv(n+l)B.(n)  +  AtD.fn)'  V(n)D.(n )]“1 

_  19.49) 

.  {Bjtn)'  [g(n+l)  +  Py(n+-l)B2(n)vUJ(n)]+  ttD^n)'  VfnjDgfnJvJn)} 

Py(n)  =  [A(n)  +  B^njK^n)]'  Pv(n+l)[A(n)  +  B1(n)Ky(n)] 

+  At[H1(n)  +  Dj(n)Kv(n)]/  V(n)[H1(n)  +  D^njK^n)] 


(9.  47) 


(9.  48) 


19.  49) 


(9.  50) 
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g(n)  =  [A(n)  +  Bx  (njKyfnJl'CgOri+i)  +  Pv(n+l)(B1(n)f(n)  +  BgCnJv^n))] 

+  AtCH^n)  +  D1(n)Kv(n)]/V(n)[D2(n)vuj(n)  +  DjWftn)]  (9*  5i) 


The  gain  Kq  (n)  is  the  solution  to  the  above  where  V(n)  is  replaced  by  Q(n). 
A  major  simplification  which  has  been  found  satisfactory  in  [46  J  and  [ll]  is 
setting 


(9.51a) 


V(n)  =  Q(n),  OsnsN 
This  is  assumed  in  the  implementation  of  the  above  equations. 

The  solution  to  the  state  estimation  problem  is 

(9.52) 

(9.53) 

where  L(n)  is  obtained  from  the  solution  of  the  forward  Riccati  equation 


x(n)  =  x^n)  +  x(n) 


x^n+l)  =  (A(n)  -  AtL(n)H2(n))  x^n)  +  AtL(n)m(n) 


P^o)  =  X(o) 


AtL(n)  =  jAfnjP^nm^n)'  +  B3(n) 


[,H^(n)PTJn)H2(n)/' 


W2(n) 


W3(n) 

At 

,-l 


L 


At 


(9.54) 


W,(n) 

P  (n+1)  =  A(n)P  (n)A(n)  +  B,(n)  — —  B„(n)' 
r\  T)  3  At  3 


(9.55) 


-  AtL(n) 


H^njP^nJH^n)'  + 


W2(n) 


At 


L(n)'  At. 


The  matrix  Pjn)  in  these  equations  is  the  covariance  matrix  of  the  estimation 
error  x(n)  given  by 


x(n)  =  x(n)  -  x  (n) 


P^(n)  =  covjx  (n)x(n)'j 


(9.56) 

(9.57) 


=  E[x(n)x(n)'} 

The  estimation  error  xta)  has  zero  mean 
E|x(n)^  =  0 


(9.58) 
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A  state  covariance  transition  matrix  can  be  derived  from  this  property. 


With 

x(n+l)  -  A(n)x(n)  +  B.(n)K(n)x(n)  +  B„(n)v  (n)  +  B„(n)ri1(n) 

1  ^  w  d  1  (9#  59) 

x(n+l)  =  A(n)x(n)  +  B1<n)K(n)x  (n)  +  BgOtfVjn) 

then 

x(n+i)  -  x(n+l)  =  A{n)(x(n)  -  X  (n))  +  B1(n)K(n)(x(n)  -  3£(n))  +  B3(n)ri1(n) 

=  (A(n)  +  B.(n)K(n))(x(n)  -  x(n))  +  A(n)x(n) 

A  (9.60) 

+  B3(n)ri1(n) 

Then 

A 

cov  {x(n+l)x(n+l)'{  =>  E  j(x(n+l)  -  x  (n+l))(x(n+l)  -  x  (n+1))'}  (9.61) 

=  (A(n)  +  BjMKfnJjE  |(x(n)  -  x(n))(x(n)  -  x(n))'}  (A(n)  +  BJ(n)K(n))/ 
+  (A(n)  +  B^nlKWJEj^n)  -  x(n))x(n)'}  A(n)' 

+  (A(n)  +  B,(n)K(n))E  |(x(n)  -  x(n))ri1(n)/}  B3(n)' 

+  A(n)E{x(n)(x(n)  -  x(n))'j  (A(n)  +  B^njKvn))' 

+  A(n)E|x(n)x(n)/}  A(n)7 
+  A(n)E{  x  (n^^n)"!  B„(n)' 

+  B3(n)E  |r^(n)(x(n)  -  x(n))'|  (A(n)  +  B1(n)K(n))/ 

+  B3(n)E  {^(njixfn)'}  A(n)' 

+  B3(n)E  |ri1(n)ri1(n)/|  B^n)' 

Since  x(n),  x(n),  and  her  e  x(n),  are  functions  of  past  inputs,  they  are 
independent  of  the  current  (white)  input  ri^n),  and 

E  |(x(n)  -  x  (n))Ti1(n)/|  =  E  jx  (n^n)'}  =  0  (9.  62) 

From  the  above 

E{(x(n)  -  x  (n))x(n/ }  =  E{x(n)x(n/  }  -  x(n)E{  x(n/ }  (9.63) 

=0-0=0 
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Hence 

cov{x(n+l)x(n+l)f  =  (A(n)  +  B^njKfn^E  |(x(n)  -  x(n))(x(n)  -  x(n))^ 

(A(n)  +  BJ(n)K(n))/  +  A(n)E|x(n)x(n)/ 1  A(n)'  +  Bg(n)E  { ^(n)  ^9‘ 64^ 
T^n/  }  B3(n) 

With  x  -  x  =  x  -  x  -  x,  (see  Figure  49. ) 

E  { (x(n)  -  x  (n))(x(n)  -  x(n)'}  =  E  |(x(n)  -  x(n))(x(n)  -  x(n))'} 

-  E  j (x(n)  -  x(n))x(n)'f  -  E{  x(n)(x(n)  -  x  (n»'f  (9.  65) 

+  Elxlnlxfn)'} 

=  cov|x(n)x(n)/}  -  2e{  x(n)x(n)/|  +  e{  x(n)x(n)/} 


cov{x(n+l)x(n+l)'}  =  [(A(n)  +  B1(n)K(n»] 

•  [cov{  x(n)x(n)/ }  -  P  (n)](A(n)  +  B.  (n)K(n))' 


+  A(n)Pt|(n)A(n)/  +  B3(n) 


Wi(n) 


(9.66) 


B3(n)  . 


Thus  the  optimal  state  covariance  matrix 


X(n)  =  E{Cx(n)  -  x(n)Ix(n)  -  x  (n)]'f 


satisfies  the  difference  equation 


X(n+i)  =  [A(n)  +  BjK(n)]  [X(n)  -  P^nJlCAfo)  +  B^n)]' 
+  A(n)PT|n)A(n)/  +  (At)"1B3(n)W1(n)B3(n/ 


(9.67) 


with  X(o)  ^  X, 


The  response  covariance  matrix,  S(n)  =  E{[r(n)  -  r(n)][r(n)  -  r(n)]  ), 
may  be  obtained  as  follows: 

r(n)  -  r(n)  =  H^nJCxfn)  -  x(n)]  +  D^njKinJCx  (n)  -  x(n)] 

=  [HjCn)  +  D1(n)K(n)]Cx(n)  -  X(n)]  -  D1(n)K(n)[x(n)  -  x(n)] 
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S(n)  =  [H^n)  +  DjInlKWlXWCH^n)  !-  D^njKtn)]' 

+  D^(n)K(n)P^(n)K  (nfr^  (n)'-  [^(n)  +  D^njKfn)] 

•  E  |[x(n)  -  X(n)][(n)  >  x(n)]'|  K'(nfD1  (n)'-  D^njKfn) 

•  E{  Cx<n)  -  x(n)][x(n)  -  x(n)yf  [H^n)  +  D^njKln)]' 

=  [H^n)  +  D^KMlXMCH^n)  +  DjMKfo)]' 

+  D1(n)K(n)PT1(n)K  (njb^n)  -  [H^n)  +  D^njK^)] 
•P^njK  (n/Dx  (n)'-  D^njK^P^^CH^n)  +  Dj(n)K(n)]/ 
Thus  the  response  covariance  matrix  is  given  by 


S(n)  =  [H^n)  +  DjInjKfnJlXfn)  -  P^nttCH^n)  +  D?(n)K(n)]' 
+  HjtaJP^nJHj  (n)/ 


(9.68) 


For  the  special  case  in  which  it  is  assumed  that  the  complete  state  can 
be  measured  exactly,  m(n)  =  x(n)  and  the  above  results  are  simplified  since 
x(n)  =  x(n)  and  P^n)  =0.  The  mean  optimal  response  vector  is  obtained  by 
substituting  (9. 44)  into  (9.36)  and  averaging  the  resulting  equation.  It  is 
given  by 


r(n)  =  H1(n)x(n)  +  Dx(n)[K(n)x  (n)  +  f(n)]  +  D2(n)va[(n)  (9.  69) 


where 


x(n+l)  =  A(n)x(n)  +  B1CK(n)x(n)  +  f(n)]  +  B2(n)viaj(n)  (9.70) 


This  finishes  the  discussion  on  the  development  of  the  optimal  control 
and  estimation  algorithms.  These  algorithms  are  implemented  as  program 
DISCOP.  The  discretized  dynamics  of  the  overall  system  is  shown  in  Figure 
50. 

The  gains  and  performance  values  obtained  from  the  equations  are 
functions  of  sample  time  At.  As  At  goes  to  zero,  the  values  obtained  by 
DISCOP  approach  to  the  continuous  model  solutions. 
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SECTION  X 


DEVELOPMENT  OF  A  METHOD  FOR  STATIONARY 
OPTIMAL  CONTROLLER  DESIGN 


As  developed  in  Section  IX,  the  computation  of  optimal  controllers  involves 
integration  of  Riccati  differential  equations  backward  for  the  controller  gains 
and  forward  for  the  estimator  »ains  together  with  state  covariance  differential 


equations. 


If  the  dynamics  are  stationary,  and  constant  gains  are  ust'i.  very  substan¬ 
tial  savings  can  be  achieved  by  directly  computing  steady-stave  '  .lu  ons  of  the 
covariance  and  Riccati  equations  [5 7} 


In  the  following,  the  development  of  stationary  design  equations  and  the 
inscription  of  algorithms  for  solving  these  equations  are  briefly  presented. 


DESCRIPTION  OF  ALGORITHM  LYAK 


The  LYAK  is  an  iterative  algorithm  for  solving  either  of  the  following 
matrix  equations  for  the  unknown  matrix  X  given  the  matrices  A  and  Q 


XA  +  A'X  +  Q  =  0 
XA'  +  AX  +  Q  =  0 


(10. 1) 
(10.2) 


In  what  follows  the  method  for  solution  is  briefly  stated.  Next  the  con¬ 
vergence  criteria  is  explained. 


Method  of  Solution 


The  method  used  to  solve  the  equations  is  iterative  and  based  on  conformal 
mapping  and  matrix  functions  [53,  57].  Given  the  matrices  A  and  Q  where  A  is 
a  stability  matrix  (real  parts  of  all  eigenvalues  of  A  are  neg^iive),  let  i//  = 
(A-arl)-l  where  or  is  a  positive  consta  .  v  i  1,  define 
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=  I  +  2a 

(10.3) 

j 

Xq  =  2a\js'GM/ 

(10.4) 

i 

i 

vlerative  algorithm  is  given  by  the  following  set  of  equations: 

•1 

AX.  =  0/X.^i 

(10.5) 

j 

XH1  =  V^i 

(10.6) 

•H 
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Choice  of  the  Parameter  a 

It  is  shown  in  [53]  that  there  is  an  optimal  value  for  the  parameter  n 
(i.  e.,  optimal  in  the  sense  that  the  algorithm  will  converge  in  a  minimum 
number  of  iterations),  call  it  a*.  Calculation  of  or*  requires  solving  for  the 
eigenvalues  of  A.  This  cannot  be  considered  because  it  is  too  expensive  com¬ 
putationally.  It  is  also  shown  in  [53],  however,  that  a  good  suboptimal  choice 
of  a,  call  it  a,  is  the  arithmetic -mean  of  the  eigenvalues  of  A.  Since  the  sum 
of  the  eigenvalues  of  A  is  just  the  trace  of  A 

a  =  |  tr  [A]  j/N  (10.8) 

where  tr  is  tl  e  trace  operator  and  N  is  the  order  of  the  matrix  A. 


Converge  nc  e  Or  iteria 


The  convergence  criteria  for  this  algorithm  is  a  ratio  test  which  is  per¬ 
formed  at  the  end  of  each  iteration.  The  absolute  value  of  the  ratio  AX^/Xi+i, 
for  each  element  in  the  upper  triangle  of  the  two  matrices  AX(  and  Xi+i,  is 
tested  to  see  if  it  is  less  than  or  equal  to  some  small  constant  e .  (The  value 
of  e  currently  being  used  is  0.  01. )  If  this  test  is  passed  Xi+j  is  accepted  as 
the  converged  solution.  If  the  test  is  not  passed  the  iterative  process  will 
continue. 


DESCRIPTION  OF  ALGORITHM  DL4K 

The  DIAK  is  a  doubly  iterative  algorithm  for  solving  the  algebraic  Riccati 
equation 

PA  +  A'P+PEP+Q  =  0  (10.9) 

where 

A  =  (A-EP)  (10.10) 

In  the  following,  first,  equivalence  relations  are  developed  between 
(10.  9)  and  the  optimization  problem  posed.  Then  the  method  of  solution  is 
given. 


Stationary  Optimization  Problem  for  Controller  Gains 

Given  the  time -invariant  matrices  F,  Gj,  G^,  H,  D,  Q  defining  the  con¬ 
trolled  system 


[  Uk  U JSi  ,-w 


'■”7 


where  r  is  the  vector  of  controlled  responses,  u  is  the  control-input  vector, 
and  riis  white  noise 

E{n1(t)n1(T),3  =  Wjfitt-r)  (10.13) 

Let  the  cost  of  control  be 

J  =  E{r  yQr }  (10.14) 

where  Q  2s  0,  D*QD  >0,  (F,  G^)  controllable  and  (F,  H)  observable  [47] . 

The  problem  is  to  find  the  gain  matrix  K  such  that  the  controller  u  =  Kx 
will  minimize  the  cost  J. 


The  covariance  equation  for  this  problem  is 


where 


0  =  (F+C  C)X  +  X(F+G1K)+ GgWjGg' 


X  =  E fax']  is  the  covariance  matrix 


(10.15) 


R  =  Ejrr'}  =  (H+DK)  X(H+DK)  'is  the  response  covariance  matrix 

(10. 16) 

The  cost  is 


J  =  tr  {(H+DK) '  Q(H+DK)X] 


(10.  17) 


where  tr  is  the  trace  operator.  Appending  the  covariance  equation  to  J  via 
the  Lagrange  multipliers  P  yields  the  Hamiltonian 


)H  ]  =  tr  ((H+DK) '  Q(H+DK)X}+  tr  {PC(P+G,K)X  +  X(F+G,K)' 
■  GjWjGj']} 

Taking  derivatives 


(10.18) 


(F+GjIOX  +  X(F+G,K) '  +  GgWjGg' 

(10.  19) 

-  ■ 

(F+G1K)/P+  P(F+GjK)+  (H+DK)'Q(H+DK) 

(10.20) 

0- 

{D'Q  (H+DK)  +  Gj'PjX 

(10.21) 
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The  optimal  controller  K,  for  this  problem,  is  the  solution  to  the  pair  of 
equations  (10.20)  and  (10.21).  Solving  (10.21)  for  K  results  in: 

K  =  -(D'QD)”1  [G1/P+  D'QH]  (10.22) 

Substituting  this  into  equation  (10.20)  yields  the  Eiccati  equation  in  P 

0=  (F-Gj[D  /QD]~1D  <^H)'P+  PfF-Gj  [D  #QD]"*1D  'QH) 

-P  {G1[D'QD]~1Gj'}P+  H'QH-H'QD  [D'QD]-1  D'QH  (10.23) 

Equating  matrices  in  equations  (10.23)  and  (10.9)  gives  the  following  rela- 
tionshios: 

A  =  (F-Gj  [D'QD]"1  D'QH)  (10.24) 

E  =  (Gj  [D'QD]-1  G^)  (10.25) 

Q  =  H'QH  -  H'QD  [D'QD]"*  D'QH  (10.26) 

Method  of  Solution 

Equation  (10.  9)  is  rewritten  as 

PA  +  A'P+  Q  =  0  (10.27) 

where 

A  =  (A-EP)  (10. 2*) 

Q  *  (Q+PEP)  (10.29) 

Starting  with  PQ  such  that  A  =  (A-EPC)  is  a  stability  matrix,  equation 
(10.  27)  is  solved  by  the  iterative  algorithm  LYAK.  The  solution  to  (10.27) 
is  substituted  into  (10.  28)_and  (10.  29)  and  then  equation  (10.27)  is  solved  again 
for  the  updated  values  of  A  and  [54}  This  process  is  continued  until  two 
successive  solutions  to  equation  (10.  27)  are  the  same  to  a  certain  number 
significant  figures  (i.  e. ,  the  same  convergence  test  as  is  used  in  the  algorithm 
LYAK. 
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STATIONARY  OPTIMIZATION  PROBLEM  FOR  THE 
ESTIMATOR  GAINS 

Consider  the  following  time -invariant  plant,  measurement  and  estimator 


equations: 

x  =  Fx  +  G^u  +  G2n1  (10.30) 

m  =  H2x  +  ri2  (10.  31) 

x  =  (F-LH2)  x+  Gxu+  Lra  (10.32) 

where  -qj  and  q2  are  stationary  white  noises  with 

Efri^t)  tj^t)')  =  W^tt-T)  (10.33) 

E{q2(t)  =  W26(t-T)  (10.34) 

E  { Ti  j(t)  ti2(t)  =  W35(t-r)  =  0  (10.35) 

and  L  is  the  estimator  gain  matrix,  (F,G2)  controllable  and  (F,Hg)  observable. 
Define  the  estimation  error  to  be 

x  ~  x  -  x  (10.  36) 

Then  the  differential  equation  of  the  estimation  error  is  obtained  from 

k  =  k  -  x  (10.37) 

Substituting  (10.  14)  and  (10.  16)  into  (10.  37)  yields 

x  =  (F-LH2)x+  G2P1  -  Lq2  (10.38) 

The  covariance  of  the  estimation  error  is  given  by 


P_  =  (F-LHJP  +  P.  (F-LH)'  +  LW.L#  +  G-W.G  '  P  (o)  =  X(o) 
q  t  t\  <5  &  i  &  q 

(10.39 

The  minimization  of  P^with  respect  to  L  yields  optimal  estimator  gain  as 


Substituting  (10.40)  into  (10.  39)  results  in 


where 


P„  =  (F-P  i)P  +  P  (F-P  I)'  +  P  iP  +  G_W1G/_ 

h  ti  r\  n  n  ti  ti  2  l  2 


i  =  h'2w;‘h2 


(10.41) 


(10.42) 


is  called  the  "information  rate". 

The  steady-state  value  of  the  estimation  error  covariance  is  given  by 


where 


AP+PA'+Q=0 
m  r\ 


A  =  (F  -  P  I) 

Q  =  (G2WlG^  +  P^) 


(10.43) 


(10.44) 


(10.45) 


Note  that  the  set  cf  equations  (10.43),  (10.44)  and  (10.45)  have  the  same 
structure  (i.  e. ,  duals)  of  the  equations  (10.27),  (10.28)  and  (10.29).  There¬ 
fore,  the  steady-state  estimation  covariance  and  the  estimator  gains  are 
obtained  also  by  using  the  algorithm  DIAK. 

THE  STEADY-STATE  COVARIANCE  WITH  THE 
OPTIMAL  ESTIMATOR 

The  covariance  of  the  controlled  system  with  the  estimator  is  obtained 
from  the  definition  given  by  (10.36 ): 


x  =  x  +  x 


(10.46) 


where 


x  =  the  state  of  the  estimator  dynamics 


x  =  the  estimation  error 


Clearly, 


where 


X  =  X  +  P  +  Y+Y' 
n 


X  =  E{xx#},  X  =  E{xx'},  P^  =  E{xx '}  andY  =  E{xx'} 


(10.47) 
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We  remark  here  that  for  the  optimal  estimator  gains 
Y  =  Y'  =  0 
so  that  (10.47)  reduces  to 
X  =  X+ 

Now  using  the  feedback  law  given  by 
u  =  -Kx 

and  substituting  (10.  31)  into  (10.  32)  yields 
x  =  (F-G1K)x+  LH2x+ Ln2 
From  (10.51)  and  (10.38)  one  obtains 


(F-G.K)  LH0 

'  °  ' '  5c, 


L^2 
>1  “  ^ ? 


(10.48) 


(10.49) 


(10.50) 


(10.51) 


(10.52) 


This  yields  the  following  set  of  differential  equations: 

X  =  (F-GjKlX  +  X(F-GjK) '  +  (LHgY  '  +  YHg'  L')  +  LW2I/ 
Y  =  (F-GjK)Y  +  Y(F-LH2)'  -  L(WgL  '  -  HgP  ) 

Pt,  =  <F-^)VPt|<F-LH2)/+LW2L'+G2WlG2/ 

The  initial  conditions  are  given  by 
X(0)  =  0 

Y(0)  =  E{x(0)  x'  (0)  }  =  0 

PJ0)  =  X(0) 
h 


(10.53) 


(10. 54) 
(10.  55) 


(10.56) 

(10.57) 

(10.58) 


On  the  account  of  the  optimality  the  last  term  in  (i0. 54)  vanishes  so  that  with 
equation  (10.57)  one  concludes  that 


Y(t)  —  0 


(10.  59) 
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Also  noting  that 


LW0L'  =  P  iP 

2  ri  n 

equation  (10.53)  can  be  written  as: 

X  =  (F-G,K)X  +  X(F-G.K)'+  P  IP 
1  1  q  n 

A 

The  steady-state  value  of  X  is  obtained  from 

^  A  A  ~  ~ 

AX  +  XA  +  Q  =  0 

where 

X  =  (F-GjK) 

Q  =  P  IP 
h  h 


(10.  60) 

(10.6?) 

(10.62) 

(10.63) 

(10.64) 


Once  P^and  X  are  found  X  is  obtained  from  (10. 49). 

An  alternate  way  of  computing  X  can  be  developed  by  writing  (,0r5  as 

*Pn  =  FPq+  PqP  '  -  VV  G2W1G2  '  <10*65> 

This  can  be  written  as 


P 

h 


=  (F-G1K)Pti+  Pti/^'-G1K)/+  GjKPt1+  PT)K/G1/ 
-  P  ip  +  G„ W.  G„  ' 

T)  T|  <:  1  £ 


(10.66) 


Summing  (10.61 )  and  (10.  66)  yields 

X=(F-G1K)X  +  X(F-G1K)/+  G0WnG  '+  [G.KP  +  P  K'G/]  (10.67) 

l  t  &  l  &  1  q  q  1 

Equation  (10.  67)  indicates  that  the  estimator  error  covariance  acts  as  a 
driver  in  the  state  covariance  differential  equation.  The  steady- state  value 
of  X  is  computed  from 

AX  +  XA'+Q  =  0  (10.88) 
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where 


A  =  (F-GjK) 


Q  =  (G^^'+G^P+PK'G^) 


(10.69) 
(10.  70) 


A  In  ADAPS,  the  total  covariance,  X,  is  computed  from  (10.49)  by  solving 
X  from  equations  (10.  62),  (10.  63  and  (10.  64)  using  the  algorithm  LYAK. 


SECTION  XI 

CONCLUSIONS  AND  RECOMMENDATIONS 


The  overall  objectives  of  this  study  were  threefold:  (1)  development 
of  theoretical  analyses  and  mathematical  models  for  precision  weapon 
delivery,  (2)  development  and  documentation  of  computer  analysis  programs, 
and  (3)  demonstration  o.'  -heir  use. 

These  objectives  were  primarily  met.  In  the  following,  qualitative  re¬ 
sults  and  recommendations  for  future  studies  pertaining  to  the  work  reported 
in  this  volume  arc  discussed. 


SIGNIFICANT  QUALITATIVE  RESULTS 

The  following  aspects  of  the  technique  developed  in  this  study  are  con¬ 
sidered  significant: 

The  stochastic  formulation  of  the  weapon  delivery  problem  is  meaningful 
and  tractable.  It  incorporates  into  the  design  the  stochastic  nature  of  the 
incident  winds,  the  time-varying  aircraft  and  weapon  dynamics,  and  the 
finite-time  nature  of  the  weapon  delivery  control  problem.  It  develops  the  full 
impact  error  covariance  matrix  using  the  overall  system  model.  It  handles 
high-order  system  descriptions,  arbitrary  sensor  arrangements,  arbitrary 
sensor  noise  levels,  and  arbitrary  control  points.  The  formulation  defines 
an  optimum  controller,  and  it  provides  a  criterion  for  measuring  the  quality 
of  any  linear  controller.  Its  basis  minimizing  the  CEP  at  impact  is  a  mean¬ 
ingful  and  appealing  design  motivation.  The  technique  makes  the  physical 
nature  of  the  weapon  delivery  problem  evident.  The  release  covariance  of  tne 
airframe  and  the  impact  propagation  matrix  of  the  weapon  show  where  the  con¬ 
trol  and  measurement  emphasis  Aould  be  placed  for  the  best  delivery. 


RECOMMENDATIONS  FOR  FUTURE  ANALYSIS  AND  MODELING  WORK 

Many  interesting  issues  came  up  in  the  course  of  the  study: 

•  Optimal  steering  of  aircraft  from  an  arbitrary  target  acquisi¬ 
tion  point  to  weapon  release.  This  can  be  posed  as  an  optimal 
control  problem  with  a  nonlinear  cost  functional.  It  can  be 
treated  using  iteratively  quadratic  control  in  the  quadratic  equi¬ 
valence  theory  of  Skelton. 

•  Although  no  simulation  was  required  in  this  study,  the  developed 
model  can  simulate  nonlinear  aircraft  and  nonlinear  weapons. 

To  increase  the  simulation  capability  to  automatic  weapon 
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delivery  a  model  of  the  nonlinear  weapon  release  state  predictor 
based  on  the  current  state  of  aircraft  should  be  developed  end 
incorporated  into  the  program. 

•  For  efficient  piloted-weapon  delivery  analysis,  the  model  of  the 
operator  should  be  incorporated  into  the  program. 

•  The  nominal  trajectories  for  the  linearization  are  generated 
by  using  a  soft  autopilot.  The  selection  of  suitable  gains  in  the 
autopilot  depends  on  past  experience  and  trial -and -err or  pro¬ 
cess.  To  increase  the  versatility  of  the  program,  the  algorithm 
designed  for  algebraic  trim  should  be  programmed,  tested  and 
incorporated  into  the  program. 

•  In  this  study,  an  exact  probability  density  function  of  the  overall 
system  is  developed  as  a  function  of  time.  Therefore,  the  uni¬ 
versally  used  assumption  of  small  correlations  in  the  CEP 
evaluation  should  be  removed  by  using  this  density  function.  A 
good  approximation  to  HPA  should  be  developed  using  the  full 
covariance  matrix. 


CONCLUSIONS 

A  reasonably  powerful  technique  for  the  analysis  and  design  of  precision 
weapon  delivery  systems  is  developed  in  this  study.  The  technique  employs 
nonlinear  modeling,  linearization,  stochastics,  quadratics  and  a  significant 
amount  of  digital  computation.  The  optimum  controller  it  produces  minimizes 
the  CEP  at  impact.  Optimal  time-varying  as  well  as  time -invariant  gains  can 
be  evaluated  for  various  airframes,  control  points,  measurement  points, 
and  weapons. 

The  value  of  the  approach  lies  in  its  mathematical  models  and  algorithms. 
They  provide  total  system  analysis  and  design  capability  by  a  digital  com¬ 
puter. 


APPENDIX  I 


DEVELOPMENT  OF  A  METHOD  FOR  FIXED-FORM 
OPTIMAL  CONTROLLER  DESIGN 


A  method  is  presented  in  this  appendix  which  designs  optimal  algebraic 
controllers  with  limited  number  of  states  (i.  e.,  fixed-form  controllers)  [55]. 

The  method  is  based  on  the  concept  of  orderly  reduction  of  the  elements 
of  optimal  full  gain  matrix  to  zero  for  the  states  which  are  not  measured.  This 
is  achieved  by  adjusting  the  remaining  elements  (i.  e. ,  elements  corresponding 
to  measured  states)  while  maintaining  the  optimality.  It  tacitly  assumes  the 
existence  of  solutions. 

The  program  which  implements  the  technique  is  called  "PROGRAM  PAPS". 
It  enables  one,  among  other  things,  to  asses  the  performance  degradations 
which  occur  when  the  complexity  of  a  full  set  of  optimal  gsins  and  a  Kalman 
filter  cannot  be  permitted. 

In  the  following,  the  problem  statement  is  given  first.  The  method  of 
solution  is  treated  next. 


PROBLEM  STATEMENT 

Given  a  time -invariant  stochastic  control  system 

x  =  Fx  +  GjU  +  ^  (I- 1 ) 

r  =  Hx  +  Du 

x  =  vector  of  state  variables 
u  =  vector  of  controls 
r  1  response  vector 
§  =  vector  of  disturbances  such  that 

E[?}=0  E{?(t)  5#(t)}  =  N5(t-T) 

The  problem  is  to  minimize  the  performance  index 

JCK^K3)  =  {tr  (JH+Dtt^+K3)]  Q  CH+DO^+K3)]'  xx'}  (1-2) 

with  a  controller  of  the  form 
=  (K1  +  K3)x 


u 


(1-3) 


In  '1-3),  K  is  an  (mxn)  matrix  with  the  following  zero  elements: 


(K  /.j  -  0  je  Q,  i  =  1,  2, . . . ,  m 


The  set  Q  denotes  a  prespecified  collection  of  iniegers  which  define  the 
unmeasured  states.  is  an  arbitrary  fixed  matrix. 


METHOD  OF  SOLUTION 


a  gain  matrix  K  corresponding  to  full-state  measurement  can  be  decom¬ 
posed  into  the  following  components: 


K  =  K1  +  \K2  +  K3,  X  =  l 


X  is  a  scalar  parameter  with 
K 


1  =  J  K.j  (ij)  eCl1 
■*  \  0  (ij)  4  Qj 


V 


l  Kij  (ij>  ef!2 

o  (ij) 


\ 


(1-4) 


K. . 
U 


r 


<  K..  (ij)  e  Qg 

jo  (ij)  4ci3 


where  the  sets  0 Q2,  and  C13  den^ce  preselected  collections  of  integers 
which  define  the  row  and  column  indices  of  Kl,  and  K^. 

The  necessary  condition  for  the  optimality  of  K*  is 

=  0  (1-5) 

X  =  0 

To  express  K*  as  a  function  of  X,  (1-5)  can  be  written  as 

r  3  %  ^  /  .,1  r^2 .1  r,3 


— ^  J(KX  +  XK2  +  K3) 
dK1 


■  y  J(K^  +  XK2  +  K3)  =  0  (with  K2  and  K3  constant  and  X  arbitrary) 
3K' 


which  implies  K*  =  K1^). 
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Then  by  the  Implicit  Function  Theorem  KA(X)  is  defined  by  the  following 
differential  equation: 


*• 

dK*(M  =  _  S2J  (K*  + 
^  L  aKxdX 


xk2  +  k3) 


52J(K1  +  XK2  +  K3) 

an1ax 


(1-6) 


A  solution  to  equation  (1-1)  can  be  obtained  by  starting  with  any  known 
terminal  condition  K  =  K*  +  Xk2+k3  for  X  =  l  and  integrating  (1-6)  backwards 
to  X  =  o.  In  the  program  the  terminal  condition  used  is  the  global  optimum 
of  the  performance  index  J  (i.  e. ,  the  solution  of  the  perfect  sensing  optimal 
quadratic  control  problem).  In  order  to  integrate  (1-6),  one  must  develop  the 
indicated  partials  and  select  a  numerical  integration  algorithm. 


NUMERICAL  INTEGRATION  ALGORITHM 

The  numerical  integration  algorithm  used  to  solve  (1-6)  is  a  predictor- 
corrector  scheme  [55]  which  employs  the  following  equations: 


(Predictor) 


K  +  24  55  dX  ~  59  *dX*  ^XK-1) 


+  37  <XK-2>  -  9  f1  (XK-3> 


(1-7) 


(Corrector) 


k!„  kP  |a2J(KP+XK+1K2+K3) 

‘  K+1>  = 


where  X  =  1;  X  =  X  +  KAX,  (K  =  1,  2,  ,  tt). 


aj(Kp+xK+lK2+K3) 


(1-8) 


o  x, -K  -o  .  — mx  x,  ^x,. 

In  order  to  obtain  enough  "back  information"  to  use  (1-7),  the  first  three 
predictor  steps  employ  the  following  equation 


KP  =  K1  (XK)  +  AX  3jjL  (Xr)  k  =  1,  2,  3 


(1-9) 


The  major  computational  task  is  to  evaluate  the  first  and  second  partial 
derivatives  in  (1-6)  and  (1-8).  A  method  to  evaluate  these  partials  is  given  in 
the  next  subsection. 
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METHOD  OF  EVALUATING  PARTIAL  DERIVATIVES 


Let  X  denote  the  covariance  matrix  of  system  (I- 1)  with  the  controller 
u  =  Kx 


(F+GK)X  +  X(F+GK)/  +  N  =  0 

and  corresponding  to  X,  define  an  adjoint  matrix  S  as  follows: 
(F+GK)  'S  +  S(F+GK)  +  (H+DK)  'Q  (H+DK)  =  0 


(MO) 


(Ml) 


The  performance  index  J(K)  of  (1-2)  is  given  by 
J(K)  =  tr  {(H+DK)'  Q(H+DK)X]  +  tr{SN] 


(M2) 


The  first  partial  of  J  [55]  is: 


1L-  =  2 tr  {[H+DK] 'QD+ SG  Ei;iX} 


(1-13) 


aKij 


where 


Eii  - 

13 


The  second  partials  of  J  are 


dK. . 


jfc  ■  4D'QD,^V  +  2  C(H+DK)'«D+SGlai  (^j 

l  "j  3a 

+  S  [(H+DK)'  QD+SG]J^-J  \ 
a  ^  '  ma  • 


(1-14) 


dX 


where  -  is  defined  by 
dKij 


(F+GK)  t|^-+  (F+GK)7  +  (GE1‘*)X  +  X(GE^) '  =  0 

ij  ij 


(1-15) 


J&L- 

SK..3X 


2TR  |(KD2)7  QDE%+ [(H+DK)/QD  +  SG]|k25^-+  Eij  )j 


(1-16) 
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APPENDIX  II 

DEVELOPMENT  OF  THE  NOMINAL  RELEASE  EQUATIONS 


In  the  analysis  and  design  of  weapon  delivery  systems,  for  a  given  attack 
maneuver  and  a  slant  range,  the  prediction  of  the  nominal  weapon  release 
time  is  needed. 

With  a  high-power  airborne  computer,  the  prediction  of  release  time  for 
hitting  a  target  can  be  based  on  the  six-degree -of- freedom  weapon  trajectory, 
computed  on  line,  taking  into  account  all  miss-producing  effects,  and  the 
current  states  of  the  aircraft.  This  may  be  referred  to  as  ''release  with  a  per¬ 
fect  computer. "  To  reduce  demand  on  high  computing  power,  however,  a 
simplified  model  is  usually  used. 

Since  no  real-time  simulation  is  involved  in  ADAPS,  the  prediction  of 
nominal  release  time  is  based  on  the  integration  of  six-degree -of- freedom 
trajectory  equations  (i.  e. ,  perfect  computer).  For  determining  the  magnitude 
of  timing  errors  in  release,  the  release  model  must  be  developed  separately 
by  the  user  as  it  depends  heavily  on  the  fire  control  system  being  evaluated. 

In  the  following  one  such  model  is  developed  for  completeness. 


DEVELOPMENT  OF  NOMINAL  WEAPON  RELEASE  EQUATIONS 

It  is  assumed  that  the  position  vector  of  the  bomb's  trajectory  is  given  by 

T  T  s 


rw(T)  " 


r  + 
r 


/  7rdT+  /  / 


v  dsdf 


(H-l) 


where 


o  o  o 

weapons  velocity  at  t  =  0 
weapons  acceleration 


The  acceleration  term  can  be  decomposed  into  contributions  due  to  gravity, 
aerodynamic  drag  per  unit  mass,  lift  per  unit  mass,  rotational  effects  per 
unit  mass,  etc.  In  order  to  avoid  line  integrals  given  by  (II- 1),  it  will  be 
assumed  that  all  these  effects  can  be  lumped  into  the  form  of  [56  ] 

ge  =  kg  (H-2) 
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d 


. 


That  is,  the  bomb  ’'sees*1  gravity  as  equivalent  to  g<>.  In  addition,  the  drag 
terms  affecting  the  forward  velocity  are  neglected.  With  these  assumptions. 


(II- 1)  can  be  written  as 

rw(6,  t)  =  rr  (0)  +  vr  (0)t  +  j  ge  t2  (II-3) 

— ♦ 

Referring  to  Figure  51,  the  miss-vector,  r,  (i.  e. ,  the  position  vector 
from  target  to  weapon)  can  be  expressed  as 

r(0,  'i)  =  r  (0,  r)-rT  (II-4) 

w  l 

=  rr<0)+  vr(0)  t+  Jge  t2  -  ?T  (H-5> 

Now  the  problem  becomes  finding  0  and  t  such  that 

?  (0,  t)  =  0  (II-6) 


Equation  (n-6)  forms  the  basic  part  of  the  so-called  fire  control  equations. 

This  algebraic  problem  can  be  solved  in  various  ways.  One  approach  is 
to  evaluate 


J**  -  min  min  |  r  (0,  t)  I  (11-7) 

0  T 

which  overrides  the  questions  of  existence  of  solutions  to  (II-6). 

In  the  following  the  existence  of  solutions  to  (H-6)  is  assumed  and  a  solu¬ 
tion  algorithm  is  developed  as  given  in  [42]  (see  Section  VII  also). 

In  matrix  notation  ( II- 6 )  is  expressed  as 
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With  ejection  velocity  vg  along  r  ,  the  velocity  term  in  II-9  becomes 


ve  v 


sin  0 


■V  Ve  J  '  cos  0 


<H-10) 


On  the  other  hand,  letting 
r  sin  01 


cos  0 


,  q  =  ^  ,  and,  qg 


(H-11) 


one  can  write 


r  =  Rr(0) 


(11-12) 


Substituting  (11-10),  (II- 11)  and  (11-12)  into  (II-9),  collecting  terms  and  dividing 
throughout  by  R  yields  the  normalized  miss -vector  equation 


f  (0,  t)  = 


( e ,  T) 


f2  (0,  t) 


=  {[I+Qt]r(0)+  :j~|  t2  -  tit  3  = 


where 


qe  q 


•q  q. 


(11-13) 


(H-I3a) 


T 

"  R 


(II- 13b) 


JM  f-R  sj  |sin’,\ 


(II- 13c) 


DEVELOPMENT  OF  THE  FIRE  CONTROL  ALGORITHM 
Let  C^^jand  define 

0  ({•,  h)  =  f  (?)  -  f  (?c)  (1-h)  =  0 

For  h  =  0  and  ?  =  one  gets  0  (?0,  o)  =  0. 


For  h  =  1,  0(?,  1)  =  f(^)  =  0.  This  implies  that  if  one  can  maintain  0(^,  h)  = 
0  by  properly  choosing  the  values  of  §  while  increasing  h  from  zero  to  one, 
the  value  of  %  at  h  =  1  becomes  the  solution  vector  to  f  (^)  =  0. 


Using  the  implicit  function  theorem: 


6h  =  0 


If  (d0/d§)is  invertible  then 

|  1lM 

3TT  ld§  I  ah 


(11-15) 


(11-16) 


f~-  H  —  t  *  t  <*o> 


(H-17) 


Therefore 


-  lit  '  <?o> 


(II- 18) 


With  the  given  hQ  and  ?0  this  differential  equation  is  integrated  up  to 
h  =  1  to  get  Now 

g 

f(»)  ■  tl  +  Qt]  r(8)  +  -rf  r2-!!  (11-19) 


*  p(5>  Mlrffi 


where 
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f  f  " 
*11  12 


f21 1  f22 


=  [(1+  Qt)  P  r(0)  1  Q  r (0)4  a-r] 


(11-20) 


0  1  0 
P  =  ,  a  = 

-10  \kg 


(11-21) 


If  Euler’s  integration  algorithm  is  used,  one  can  write  from  (11-18) 


(ek'  V 


(n-22) 


where 


b  =  f  Ah 
o 


and  (11-20)  is  iterated  from  k  =  0  to  k  =  ^  =  N. 


To  establish  startup  values  for  the  fire  control  algorithm,  select 
V_ 


0  =  tan 

o 


(1-23) 


and  compute  tq  so  that  zq  =  f2(0o>  =  i$ves 


^  /2(zt  -  R) 

ro  *V - leg 


(11-24) 


Substituting  this  into  fj  in  (II-9)  gives  the  normalized  range  error  at  the 
beginning  of  the  iteration 


Xo  =  fl  (V  To)  =  Vr  To  "  XT 


(11-25) 


Thus,  following  equations  are  used  to  start  up  the  iterations: 


*  i 

t 


ci 


'^-  >  V  -Tpp^rr-  tp— >-  - - - 


.-1  Ve 


1(2  -  R) 

T 

“Kg 


and  f 


V  T 
r  o 


(11-26) 


After  hr-hsg  found  £p  so  that  (11-13)  is  satisfied,  the  predicted  time  during 
"pullup"  is  given  oy 


y+  er 


(H-27) 


Due  to  simplifications  used  in  the  modeling  of  the  weapon  trajectory,  the 
predicted  nominal  release  time  tpu  must  be  corrected.  This  correction  is 
in  the  form  of 


t  =  k  t  +  At  -  At 
r  c  pu  c 


(11-28) 


where 


k  -  algorithmic  error  correction  multiplier 

v« 

At  =  algorithmic  error  correction  bias 

V 

At  =  known  delay  in  the  release  mechanism 

Obviously,  the  actual  nominal  release  will  occur  at 

t  =  k  t  +  At 
ra  c  pu  c 


(n-29) 


The  correction  multiplier  kc  and  bias  Atc  are  obtained  from  numerical  experi¬ 
ments  with  actual  six-degree-of-freedom  weapon  model. 


■.  'If  i  W ' 
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